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Abstract 

Orbit  disposal  and  maintenance  of  aging  satellites  has  become  a  significant  concern  over  the 
past  few  years,  as  the  increasing  number  of  orbiting  objects  in  high-value  orbits  (e.g.  geosyn¬ 
chronous)  threatens  to  limit  the  launching  of  future  satellites  and  spacecraft.  Many  of  the  satellites 
currently  in  orbit,  however,  were  not  built  with  disposal  considerations.  The  DSCS  II  series,  for 
example,  was  launched  into  orbit  beginning  in  the  1970s,  and  many  satellites  are  now  without  the 
fuel  required  to  conventionally  transition  to  a  sanctioned  disposal  orbit.  In  geosynchronous  orbit, 
however,  the  largest  non-gravitational  perturbation  is  solar  radiation  pressure  (SRP).  By  adjusting 
the  position  of  the  satellite  with  a  controller  to  maximize  the  perturbing  acceleration  due  to  the 
force  of  SRP,  the  satellite  can  be  slowly  raised  into  an  appropriate  disposal  orbit.  The  results  from 
this  study,  along  with  validation  results  propagated  with  Satellite  Tool  Kit,  are  presented.  After 
making  several  simplifying  assumptions  (primarily  a  cylindrical  Earth  shadow,  flat  plate  geometry, 
and  constant  cross-sectional  area  of  the  satellite),  the  time  required  to  raise  the  modelled  DSCS  II 
F-13  satellite  400  km  into  a  disposal  orbit  is  approximately  33  years.  This  time-to-disposal  can  be 
reduced  by  using  a  larger  area-to-mass  ratio  and  more  reflective  surface  materials. 
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The  Effects  of  Using  Solar  Radiation  Pressure 
to  Alleviate  Fuel  Requirements  for  Orbit  Changing 
and  Maintenance  of  the  DSCS  II  F-13  Satellite 

I.  Introduction 

1.1  Motivation:  Using  Solar  Radiation  Pressure  as  a  Controller 

Detachment  12  of  the  Space  and  Missile  Systems  Center  (SMC  Det  12)  serves  as  the  pri¬ 
mary  provider  of  launch  capability,  spaceflight  and  on-orbit  operations  for  Department  of  Defense 
users  [27].  SMC  Det  12,  and  the  Air  Force  Space  Command  in  general,  has  continuing  interest  in 
the  area  of  orbit  disposal  and  maintenance  of  aging  satellites.  In  particular,  the  study  of  moving 
an  existing  satellite  (that  has  expended  all  of  its  onboard  fuel  supply)  into  a  disposal  orbit  is  of  in¬ 
terest,  as  is,  conversely,  de-orbiting  a  satellite  using  low-fuel  maneuvers  and  potential  lunar  flybys. 
The  first  topic  will  be  addressed  in  this  thesis;  the  second  is  left  as  future  work. 

For  high  altitude  orbits,  solar  radiation  pressure  becomes  the  dominant  perturbative  force, 
especially  for  satellites  with  large  surface  areas.  By  intentionally  aligning  a  satellite’s  orientation 
with  respect  to  the  sun,  this  perturbative  force  can  be  used  to  maneuver  the  vehicle.  This  could 
potentially  be  used  as  a  means  of  moving  satellite  into  disposal  orbits  when  there  is  insufficient  fuel 
onboard.  This  topic  has  been  of  particular  interest  to  the  Air  Force  and  Department  of  Defense 
of  late  due  to  the  increasing  number  of  in-orbit  vehicles.  As  a  result,  research  into  the  disposal  of 
satellites  (and  the  related  topic  of  collision  avoidance)  has  been  encouraged. 

An  additional  benefit,  by  having  this  technique  available,  is  that  the  mission  life  of  many  Air 
Force  satellites  might  be  extended  because  less  fuel  would  be  needed  for  station  keeping  and/or 
orbit  maintenance.  This  thesis  uses  a  representative  Air  Force  spacecraft,  a  Defense  Satellite 
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Communications  System  (DSCS)  F-13  satellite,  as  a  test  vehicle  for  demonstrating  the  ability  of 
solar  radiation  pressure  to  alleviate  fuel  requirements  for  orbit  changing  and/or  orbit  maintenance. 

1.2  Background 

1.2.1  Orbital  Perturbations  and  Solar  Radiation  Pressure.  Every  object  in  space,  whether 
in  orbit  about  a  planet  (e.g.  Earth)  or  flying  through  interstellar  regions,  experiences  small  changes 
in  its  two-body  motion  caused  by  external  forces.  These  disturbances  are  called  orbital  perturba¬ 
tions,  and  fall  into  two  main  categories:  gravitational  and  non-gravitational.  Gravitational  per¬ 
turbative  forces  are  those  involving  the  gravitational  attraction  of  secondary  bodies  (i.e.  Earth, 
Moon  or  Sun)  on  the  satellite.  Non-gravitational  perturbative  forces  are  all  the  other  types  of 
disturbances,  including  atmospheric  drag,  solar  radiation  pressure  and  Earth  albedo.  For  Earth¬ 
orbiting  satellites,  such  as  the  DSCS  II  F-13,  the  dominant  perturbations  are  gravitational.  The 
common  spacecraft  perturbations  for  satellites  orbiting  in  the  geosynchronous  belt  are  listed  in 
Table  1.1  [10]. 

In  this  thesis,  the  primary  gravitational  perturbations  that  will  be  considered  are  from  the 
Sun  and  the  Earth.  The  uneven  mass  distribution  of  the  Earth  causes  zonal  harmonics,  which  act 
to  alter  the  apparent  forces  of  the  Earth  on  the  satellite  depending  on  where  the  satellite  is  orbiting 
about  the  Earth.  The  main  gravitational  perturbation  for  a  satellite  in  geosynchronous  orbit  is 
the  second-order  even  zonal  harmonic,  commonly  referred  to  as  the  Earth  oblateness  and/or  J2 
perturbation.  The  solar  third-body  perturbation  affects  the  satellite  orbit  slightly  less,  although 
the  modelling  of  the  disturbance  is  significantly  more  difficult  computationally  than  for  the  J2 
perturbation  due  to  the  complicated  orbit  of  the  Earth  about  the  Sun  and  the  magnitude  of  the 
vectors  involved.  The  orbital  perturbations  caused  by  the  Sun’s  gravity  are  less  than  for  the  Moon 
(which  will  be  neglected  due  to  computational  instability,  with  no  diminished  value  of  the  thesis), 
but  have  a  greater  effect  than  solar  radiation  pressure. 
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Table  1.1:  Common  Geosynchronous  Satellite  Perturbations, 

Adapted  from  [10]. 


Perturbation 

Formula 

Acceleration  (m/s2)  of  Geosynchronous 
Spacecraft  with  A/M  =  0.01  m2 /kg 

Earth’s  Oblateness 

3GMe  (^e)2  J2Q 

7.4  x  10"6 

Lunar  Third  Body 

2®-r 

7.3  x  10-6 

Solar  Third  Body 

2gmq 

4 

3.3  x  10-6 

Atmospheric  Drag 

iCdwPV2 

negligible 

SRP 

A  $© 

M  c 

4.6  x  10"8 

Earth’s  Albedo 

4.2  x  lO"10 

The  primary  non-gravitational  perturbation  for  a  geosynchronous  satellite  like  the  DSCS  II  F- 
13  comes  from  solar  radiation  pressure  (SRP).  The  SRP  disturbance  is  caused  by  incoming  incident 
solar  photons  impinging  on  an  energy  absorbing  surface  (the  satellite) ,  which  generates  a  force  per 
unit  area  on  the  surface.  If  the  vector  sum  of  all  the  forces  acting  on  the  distinct  segments  of  the 
satellite  does  not  pass  through  the  satellite’s  center  of  mass,  a  torque  will  also  be  generated  [8]. 
The  SRP  perturbation  is  essentially  equivalent  to  adding  an  incremental  velocity  (Aw)  to  one  end 
of  the  orbit  (e.g.  perigee)  and  subtracting  the  same  Av  at  the  other  end  (e.g.  apogee),  causing  the 
orbit  to  gradually  become  elliptic.  Six  months  later,  the  Earth  is  at  the  opposite  side  of  its  orbit 
around  the  Sun,  and  the  action  is  reversed  to  recircularize  the  orbit  [9].  Figure  3.7  illustrated  this 
concept. 

1.2.2  DSCS  II  F-13.  The  satellite  chosen  for  study  in  this  thesis  is  the  DSCS  II  F- 
13,  which  was  launched  into  a  geosynchronous  orbit  on  21  November  1979  from  Cape  Canaveral, 
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Florida.  The  16  TRW-built  DSCS  II  satellites  form  a  constellation  of  communications  satellites 


once  charged  with  relaying  secure  voice  and  data  communications  for  the  United  States  military. 
They  replaced  a  set  of  26  spacecraft  launched  under  the  Department  of  Defense’s  Initial  Defense 
Communications  Satellite  Program  beginning  in  1966.  Each  DSCS  II  spacecraft  is  a  cylindrical 
spin-stabilized  satellite  with  a  despun  antenna  platform  housing  two  parabolic  antennas.  The  spin- 
stabilized  cylinder  is  covered  with  a  layer  of  body-mounted  solar  arrays  to  provide  electrical  power 
augmented  by  a  system  of  three  Nickel  Cadmium  batteries. 

The  DSCS  II  F-13’s  initial  mass  was  611  kg;  it  had  an  initial  perigee  height  of  37104  km  and 
an  initial  apogee  height  of  37195  km.  Upon  reaching  orbit,  the  satellite’s  inclination  was  13.6°. 
Its  design  life  was  only  intended  to  be  five  years,  but  the  satellite  provided  many  more  years  of 
valuable  service  to  the  Department  of  Defense  (DoD).  Additional  parameters  for  the  DSCS  II  F-13 
can  be  found  in  Appendix  A  [23,37]. 

1.3  Research  Goals 

1.3.1  Problem  Statement.  The  primary  objective  of  this  thesis  is  to  determine  how  long 
it  will  take  to  raise  the  orbit  of  the  DSCS  II  F-13  satellite  into  a  disposal  orbit  using  SRP.  To  meet 
this  goal,  this  thesis  undertakes  the  development  of  a  numerical  model  in  Matlab®  to  simulate  the 
satellite’s  orbit,  orientation  with  respect  to  the  sun  and  the  effects  of  the  perturbative  forces  caused 
by  the  solar  radiation  pressure,  Earth  oblateness  and  the  solar  third-body  effects.  In  addition  to 
the  simulation  effort,  the  thesis  involves  the  development  of  two  techniques  to  control  the  attitude 
of  the  satellite  by  taking  advantage  of  the  applied  solar  radiation  pressure  forces. 

1.3.2  Scope  and  Assumptions.  The  DSCS  II  F-13  satellite  will  be  modelled  with  approxi¬ 
mated  despun  and  spin-stabilized  platform  areas.  The  code  will  implement  a  case  structure,  which 
will  enable  several  parameters  to  be  easily  varied:  despun  platform  area,  coefficients  of  absorption 
and  specular  and  diffuse  reflection  of  the  despun  platform  area,  number  of  orbits  and  the  controller 
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used.  The  solution  will  be  propagated  with  a  cylindrical  Earth  shadow,  constant  cross-sectional 
area  of  the  satellite  and  the  perturbations  addressed  in  Section  1.3.1.  The  modified  equinoctial 
equations  of  motion  of  the  satellite  will  be  integrated  to  obtain  the  solution  and  the  results  will  be 
plotted. 

For  the  purposes  of  this  thesis,  the  collision  avoidance  problem  will  not  be  addressed.  Nor 
will  the  number  of  man-hours  required  to  raise  the  orbit  have  a  consideration  on  the  complexity  of 
the  controller  used,  although  it  would  indeed  in  the  operational  world.  The  contributions  of  several 
perturbations  will  not  be  examined,  due  to  the  small  effects  they  have  on  satellites  in  geosyn¬ 
chronous  orbit.  These  include  perturbing  forces  due  to  atmospheric  drag,  the  Earth’s  magnetic 
field,  Earth  albedo  and  infrared  radiation  pressure,  higher  order  even  zonal  harmonics,  sectorial 
harmonics,  tesseral  harmonics  and  tides.  The  model  will  not  include  periods  of  dramatically  in¬ 
creased  or  decreased  solar  activity  or  a  conical  Earth  shadow  model  due  to  the  small  percentage 
gain  in  accuracy  (approximately  5%  for  each) .  Since  the  unknown  surface  properties  of  the  satellite 
cause  the  results  to  vary  by  as  much  as  80%,  the  additional  computation  time  required  to  model 
these  perturbations  yields  no  real  added  value. 

To  limit  the  run  time  of  the  code  to  a  reasonable  threshold,  a  number  of  simplifying  assump¬ 
tions  were  used.  As  previously  mentioned,  a  cylindrical  shadow  model  is  used  that  neglects  the 
penumbra  and  umbra  regions.  This  shadow  model  essentially  regards  the  Sun  as  a  point  source,  and 
neglects  the  transition  regions  of  the  shadow  projection.  Additionally,  the  shadow  is  considered  to 
be  unaltered  due  to  any  flattening  of  Earth’s  pole  and/or  changes  in  atmospheric  density,  although 
both  could  be  considered  for  future  work. 

Further  assumptions  include  a  constant  (original)  satellite  mass  and  approximated  satellite 
surface  areas.  The  satellite  parameters  are  not  well  documented;  the  area  of  the  parabolic  antennae 
on  the  despun  platform  can  only  be  estimated.  For  the  purposes  of  this  thesis,  the  parabolic 
antennas  will  be  combined  and  modelled  as  one  flat  plate  with  zero  thickness.  The  mass  of  the 
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satellite  is  considered  to  be  negligible  for  purposes  of  the  solar  third-body  perturbation  calculations. 
Additionally,  the  nominal  coefficients  of  specular  reflection,  diffuse  reflection  and  absorption  are 
taken  to  be  one-third  each,  since  the  current  reflective  properties  of  the  satellite  are  not  well  known. 

The  SRP  model  used  in  the  development  of  the  simulation  code  ignores  any  torque  caused 
by  SRP  resultant  force  acting  outside  the  satellite’s  center  of  mass.  Several  algorithms  from  Val- 
lado’s  text  were  referenced  in  the  development  of  the  code,  and  they  each  have  their  own  inherent 
assumptions.  These  assumptions  are  minor  as  compared  to  the  simplifying  assumptions  already 
addressed,  and  may  be  referenced  in  the  text  accompanying  the  algorithms  in  his  book  [32]. 

For  the  purposes  of  this  thesis,  the  desired  disposal  orbit  is  defined  as  400  fcm  above  the  epoch 
altitude  of  the  satellite  for  the  DSCS  II,  as  it  already  resides  in  the  disposal  orbit  set  forth  by  the 
United  States  Space  Command  in  Appendix  D.  Additionally,  the  orbital  period  is  considered  to  be 
constant  over  the  integration  time,  and  the  equatorial  radius  of  the  Earth  is  considered  to  remain 
at  a  constant  6378  km. 

Finally,  Coverstone’s  controller  model  (see  Section  3.4.2)  has  inherent  simplifications  and 
assumptions.  Most  notably,  she  assumes  only  specular  reflection  and  a  constant  solar  flux.  Further 
information  regarding  her  assumptions  can  be  found  in  her  paper  [11]. 

1-4  Overview  of  Chapters 

The  structure  of  the  thesis  is  such  as  to  guide  the  reader  through  the  thought  processes 
involved.  The  next  chapter  summarizes  the  previous  research  conducted  relating  to  the  subject 
matter,  and  provides  the  groundwork  for  cultivating  an  appropriate  model  for  the  orbit  of  the 
DSCS  II  F-13.  Chapter  III  presents,  in  detail,  the  theory  and  equations  behind  the  algorithms 
used  to  create  the  model  of  the  orbit.  Chapter  IV  will  enumerate  the  numerical  results  obtained 
from  the  Matlab®  model  developed  in  Chapter  III.  Recommendations  and  conclusions  will  then 
be  outlined  in  Chapter  V,  along  with  suggestions  for  any  future  work  that  might  be  performed. 
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II.  Background 


2.1  The  Infancy  of  Solar  Radiation  Pressure  Research 

The  earliest  studies  of  SRP  began  in  the  late  19th  century.  Not  until  the  first  spacecraft 
were  launched  in  the  middle  of  the  20th  century,  however,  did  in  situ  data  become  available,  and 
hence  theories  were  experimentally  proven.  Over  the  past  century,  our  knowledge  of  SRP  and  its 
perturbing  effects  on  spacecraft  in  Earth  orbit  and  interplanetary  space  has  matured  and  developed 
into  an  intricate  science.  For  many  years,  the  focus  of  SRP  studies  was  confined  to  the  minimization 
of  its  effects  on  spacecraft.  More  recently,  considerable  effort  has  been  put  into  the  development  of 
applications  to  harness  the  SRP  effects  for  purposes  of  interplanetary  propulsion,  namely  with  solar 
sails.  Consequently,  the  work  begun  by  many  pioneers  of  the  late  1800s  and  early  1900s  has  come 
full  circle  as  their  theories  are  now  used  to  harness  the  SRP  acceleration  as  opposed  to  minimizing 
its  effects.  The  focus  of  this  thesis  rests  in  this  area:  using  the  perturbing  SRP  acceleration  to  raise 
the  orbit  of  a  satellite. 

Although  this  thesis  concentrates  on  the  effects  of  SRP,  it  is  important  to  note  that  the  held 
of  study  grew  from  work  regarding  radiation  pressure  on  Earth.  In  1873,  Scottish  physicist  James 
Clerk  Maxwell  set  forth  the  first  suppositions  about  the  existence  of  radiation  pressure  [29].  For 
nearly  30  years  after  Maxwell,  the  subject  matter  laid  effectively  dormant,  until  Nichols  and  Hull 
from  Darthmouth  College  and  Russian  physicist  Peter  Lebedev  were  able  to  (independently)  exper¬ 
imentally  measure  radiation  pressure  in  precision  laboratory  tests  in  1900.  About  the  same  time, 
one  of  Maxwell’s  students,  English  physicist  and  mathematician  John  Henry  Poynting,  suggested 
the  existence  of  an  effect  of  the  Sun’s  radiation  that  causes  small  particles  orbiting  the  Sun  to 
gradually  approach  it  and  eventually  enter  the  Sun’s  atmosphere. 

The  earliest  study  of  SRP  in  particular  began  in  1924,  when  Russian  astronautics  pioneer  Kon¬ 
stantin  Eduardovitch  Tsiolkovsky,  widely  considered  the  father  of  astronautics  and  rocket  dynamics, 
and  his  colleague,  Russian  engineer  Fridrickh  Arturovich  Tsander,  first  introduced  a  proposal  to 
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use  SRP  as  a  method  of  spacecraft  propulsion.  Their  concept  used  very  thin  sheets  of  mirrors  with 
large  surface  areas  to  harness  the  Sun’s  radiation  pressure  as  a  form  of  propulsion.  These  concepts 
matured  into  the  evolution  of  practical  solar  sailing,  first  written  of  by  Tsander  around  1924  [22]. 
During  this  period,  Poynting  completed  the  first  study  of  the  Sun’s  radiation  pressure  as  it  affects 
small  meteorites  in  interplanetary  space.  Later,  American  cosmologist  Howard  Percy  Robertson 
expanded  on  Poynting’s  research  and  the  resulting  theory  is  now  called  the  Poynting-Robertson 
effect. 

2.2  The  Space  Age 

After  the  first  spacecraft  (Sputnik  I  and  II,  Explorer-I  and  Vanguard)  were  successfully 
launched  in  the  late  1950s,  in  situ  data  on  the  perturbing  effects  of  SRP  became  available.  The 
launch  of  the  Vanguard  and  EcholA  satellites,  and  the  ensuing  perturbations,  fueled  scientists 
around  the  world  to  study  the  effects  closely  and  develop  accurate  models  that  would  predict  the 
effects  of  SRP  on  spacecraft  during  their  orbital  lifetime.  Koskela,  Kozai  and  Musen,  three  scien¬ 
tists  that  studied  these  early  perturbations,  each  sought  to  build  models  to  describe  the  SRP  effects 
on  satellites  in  Earth  orbit.  These  men  all  used  simplifying  assumptions  in  their  derivations,  most 
notably  constant  solar  flux,  constant  satellite  area  (perpendicular  to  the  solar  radiation  vector), 
cylindrical  Earth  shadow,  no  diffuse  reflection  and  no  absorption  [10]. 

In  the  early  1960s,  scientists  began  to  hone  the  models  developed  by  the  first  researchers 
(namely  Koskela,  Kozai  and  Musen).  Robert  Bryant,  who  was  working  for  the  National  Aeronautics 
and  Space  Administration  (NASA)  Goddard  Space  Flight  Center,  sought  to  refine  the  previous  work 
by  conducting  his  own  study  in  1961.  He  was  the  first  to  observe  that  when  Earth’s  shadow  is 
neglected  in  the  orbital  propagation  model,  a  spacecraft  will  only  experience  short  periodic  changes 
in  its  semimajor  axis.  Conversely,  when  the  SRP  model  includes  the  shadow,  the  satellite  will 
experience  a  more  significant  perturbation.  This  had  great  implications  for  users  at  the  time, 


since  the  principal  finding  of  Bryant’s  work  indicated  that  the  long-term  effects  of  SRP  are  much 
greater  when  a  satellite  experiences  solar  eclipsing  by  the  Earth  [7, 10].  This  discovery  will  prove 
to  be  a  driving  factor  in  the  development  of  this  thesis,  as  the  problem  requires  not  the  short-term 
perturbations  of  the  semimajor  axis,  but  the  long-term  effects  obtained  by  spending  time  in  Earth’s 
shadow. 

2.3  Modern  Research  and  Applications 

As  the  years  elapsed  from  Bryant’s  groundbreaking  discovery  in  the  early  1960s,  spacecraft 
rapidly  evolved  into  more  complex  vehicles.  The  equations  of  motion  that  were  sufficient  to  properly 
track  orb-shaped  satellites  in  low-earth  orbits  were  no  longer  adequate.  More  powerful  launch 
vehicles  enabled  satellites  to  reach  higher  orbits,  and  the  effects  of  SRP  began  to  dominate  the 
perturbations.  New  generations  of  scientists  were  directed  to  study  the  effects  of  SRP  under  these 
new  conditions. 

Harwood  et  al.  began  with  the  simplest  of  satellites,  a  spherical  spacecraft.  This  spacecraft 
had  a  constant  cross-sectional  area  and  an  acceleration  vector  that  was  perpendicular  to  the  SRP 
vector,  which  greatly  reduced  (with  the  assistance  of  other  simplifying  assumptions)  the  computa¬ 
tional  time  need  to  analyze  the  SRP  perturbations.  Harwood  introduced  a  varying  Earth-Sun  (and 
therefore  satellite-Sun)  distance,  which  enabled  the  solar  flux  constant  to  be  realistically  varied 
according  to  time  of  year.  For  the  first  time,  a  time-varying  SRP  acceleration  was  studied. 

The  first  satellite  to  require  precise  SRP  modeling  was  the  oceanographic  satellite  Ocean 
Topography  Experiment  (TOPEX)/Poseidon  .  This  satellite  was  built  to  acquire  altimeter  mea¬ 
surements,  which  are  used  to  map  the  ocean  topography.  TOPEX/Poseidon  justified  further  study 
and  application  of  more  accurate  models  for  SRP  disturbing  acceleration,  because  even  small  errors 
in  the  orbit  determination  process  would  be  magnified  into  significant  inconsistencies  in  its  topo¬ 
graphical  measurements.  Marshall  et  al.  developed  a  more  rigorous  modeling  tool  that  simulated 
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the  satellite  with  a  box-wing  shape  (a  box  as  the  main  body  and  flat  plates  for  the  solar  panel 
arrays).  The  satellite’s  cross-sectional  area  was  allowed  to  shift  as  it  rotated  around  the  Earth  and 
underwent  programmed  attitude  changes.  Like  earlier  studies,  some  simplifying  assumptions  were 
made  by  the  authors,  such  as  a  cylindrical  Earth  shadow  and  constant  surface  reflection  properties. 
The  authors  implement  Lambert’s  cosine  law  for  diffuse  reflection  (as  does  Vallado  [32]),  which, 
although  a  simplification  of  the  full  diffusion  model,  is  more  than  most  authors  use  in  their  research. 
In  Marshall’s  model,  the  authors  also  accounted  for  the  lesser  perturbations  due  to  the  satellite’s 
thermal  emissions,  Earth’s  albedo,  and  Earth’s  infared  emissions  [10]. 

Much  of  the  contemporary  modeling  research  efforts  are  focused  on  the  effects  of  SRP  per¬ 
turbations  on  interplanetary  spacecraft.  One  such  project  is  the  Japanese  Selenological  and  En¬ 
gineering  Explorer  (SELENE),  which  will  use  an  octagon-shaped  relay  satellite  to  transmit  com¬ 
munications  between  SELENE’s  main  orbiter  and  Earth  following  SELENE’s  launch  in  2007  (see 
Figure  2.1).  The  long-term  effects  of  SRP  on  the  relay  satellite  will  be  studied  for  approximately 
one  year,  once  it  achieves  its  elliptical  orbit  around  the  Moon  [36].  Since  the  Moons  mass  is  rela¬ 
tively  insignificant,  it  cannot  hold  much  of  an  atmosphere,  and  hence  SRP  becomes  the  dominant 
non-gravitational  perturbation  acting  on  the  satellite.  For  this  reason,  the  effects  of  SRP  on  the 
lunar-orbiting  satellite  must  be  modeled  accurately.  Project  engineers  used  force  equation  models 
similar  to  those  postulated  by  Chobotov,  Ries  et  al.  and  Milani  et  al.  in  their  derivation  of  the 
SRP  effects,  but  have  neglected  the  shadow  of  the  Moon  [10]. 

A  program  widely  used  in  the  astrodynamics  community  to  study  the  general  effects  of  per¬ 
turbations  was  developed  by  the  Charles  Stark  Draper  Laboratory.  The  Draper  Semi-analytical 
Satellite  Theory  (DSST)  model  analytically  separates  the  secular  and  long-periodic  motions  from 
the  short-periodic  motion,  thus  enabling  the  user  to  modify  the  modeled  forces  to  suit  their  accu¬ 
racy  and  computation  time  requirements.  DSST  makes  several  assumptions;  the  SRP  model  uses 
a  cylindrical  Earth  shadow  and  a  constant  coefficient  of  reflection.  The  DSST  model  is  convenient 
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Figure  2.1:  SELENE  Orbiter.  The  spacecraft  has  three  components:  the  main  orbiter,  a  small 

relay  satellite,  and  a  small  Very  Long  Baseline  Interferometry  (VLBI)  satellite.  The  main  orbiter, 
which  has  a  rectangular  box  shape,  carrys  the  scientific  instrumentation.  The  relay  satellite  is 
an  octagonal  shape  and  will  be  used  to  communicate  between  the  orbiter  and  Earth.  The  VLBI 
satellite  is  the  same  shape  as  the  relay  satellite;  its  purpose  is  to  conduct  investigations  on  the 
position  and  precession  of  Earth’s  moon  [24] . 


11 


for  many  users  because  it  tailors  the  program  to  take  advantage  of  the  best  of  both  general  and 
special  perturbation  techniques.  Cefola  et  al.  have  recently  used  the  DSST  program  to  study  the 
short-periodic  motion  of  formation  flying  satellites  [25] . 

An  additional  study  of  non-gravitational  perturbations  was  presented  by  Bruce  Bowman  et 
ah,  of  the  Space  Warfare  Center  at  Schriever  Air  Force  Base  (AFB)  in  Colorado,  in  2000.  His  paper 
addressed  the  long-term  effects  of  SRP,  atmospheric  drag  and  Earth  albedo  on  the  approximately 
20  West  Ford  needles  clusters  that  formed  after  the  initial  package  was  launched  in  1963.  The 
copper  needles,  which  have  a  large  area  to  mass  ratio,  decayed  rapidly  once  placed  in  orbit,  and 
clusters  were  formed  due  to  delays  in  releasing  the  needles. 

Bowman’s  analysis  is  of  interest,  as  the  long-term  perturbations  due  to  direct  SRP  found 
in  the  West  Ford  needles  could  be  analogous  to  the  attempted  raising  of  an  orbit  in  this  thesis. 
Although  the  needles  are  primarily  located  in  a  low  earth  orbit  (LEO)  orbit,  without  active  control 
the  West  Ford  needles  exhibited  a  semimajor  axis  displacement  of  approximately  10  km.  in  34 
years  [10].  In  a  GEO  orbit,  as  for  this  thesis,  the  semimajor  axis  should  be  expected  to  raise  much 
more  in  the  same  amount  of  time  (with  a  controller),  so  the  results  from  Bowman  are  encouraging. 

2.3.1  Earth  Shadow  Models.  It  is  important  to  note  at  this  point  that  many  researchers 
over  the  past  century  have  made  simplifying  assumptions  regarding  the  Earth’s  shadow,  namely 
by  either  using  a  cylindrical  model  for  analysis  or  neglecting  the  shadow  entirely.  For  calculations 
where  the  effects  of  SRP  on  a  spacecraft  were  not  required  to  be  known  absolutely,  or  over  long 
periods  of  time,  neglecting  the  shadow  is  acceptable.  In  the  application  addressed  in  this  thesis, 
however,  SRP  effects  play  a  major  role  in  the  outcome,  and  the  simplifying  assumptions  made  by 
previous  authors  will  not  provide  an  accurate  enough  model. 

One  of  the  simplest  models  used  to  model  SRP  in  the  presence  of  Earth’s  shadow  was  proposed 
in  1971  by  Sylvio  Ferraz  Mello.  He  introduced  a  shadow  function,  ip,  which  is  set  equal  to  one  if 
the  satellite  is  illuminated  and  to  zero  when  the  satellite  is  in  shadow  [12].  A  more  sophisticated 
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method  was  developed  by  Kabelac  in  1988.  He  derived  a  method  which  uses  a  more  realistic  conical 
Earth  shadow  with  both  penumbra  and  umbra  regions.  Kabelac’s  model  accounts  for  the  mutual 
position  of  the  Sun,  Earth  and  satellite,  the  shape  of  conical  surface,  light  diffusion,  atmospheric 
absorption  and  the  effect  of  ozone  [16]. 

A  number  of  researchers  have  used  the  updated  shadow  models  to  analyze  spacecraft  eclipse 
transition  effects  in  further  detail.  Vokrouhlicky  et  al.  have  written  a  sequence  of  papers  on  this 
topic.  The  papers  use  Kabelac’s  model  to  analyze  the  shadow  effects  on  a  spacecraft  as  it  transitions 
between  regions  of  the  shadow  projection.  The  authors  note  that  the  apparent  size  of  the  sun,  with 
respect  to  the  spacecraft,  changes  as  the  spacecraft  moves  in  and  out  of  the  penumbra  region  [33]. 
The  authors  then,  in  later  papers,  expound  upon  Kabelac’s  model  to  include  the  altering  of  Earth’s 
shadow  due  to  the  flattening  of  Earth’s  pole  and  changes  in  atmospheric  density  [10].  This  particular 
inclusion  is  beyond  the  scope  of  this  thesis  but  could  be  addressed  in  follow-on  studies. 

2.3.2  Advanced  Solar  Radiation  Pressure  Force  Models.  Modern  scientists  have  also 
focused  on  developing  more  accurate  solar  radiation  pressure  force  and  torque  models.  R.M. 
Georgevic,  perhaps  the  first  to  consider  modeling  the  effects  of  SRP  on  reflecting  surfaces,  de¬ 
veloped  a  mathematical  archetype  for  computing  the  radiation  pressure  force  and  torques  on  both 
a  parabolic  high  gain  antenna  and  a  circular  cylinder  in  1973.  Georgevic  emphasized  the  impor¬ 
tance  of  approximating  a  body’s  irradiated  surface  area  and  determining  the  specular  and  diffuse 
reflection  vectors  [13].  David  Vallado  based  his  model  for  SRP  acceleration  on  Georgevic’s  early 
work  and  on  Burns  et  al.  (2000).  Burns  et  al.  developed  an  attitude  guidance  law  for  time- varying 
spacecraft  orientation,  which  is  more  realistic  than  assuming  the  satellite  maintains  a  constant  at¬ 
titude  normal  to  the  Sun.  Vallado’s  model  includes  absorption  and  specular  and  diffuse  reflection. 
To  model  the  diffuse  reflection,  he  assumes  a  Lambertian  reflective  surface;  this  means  that  light 
impinging  on  the  surface  is  scatted  such  that  the  surface  luminance  is  the  same  regardless  of  the 
observer’s  viewing  angle.  Not  all  surfaces  are  perfect  Lambertian  reflectors,  but  the  approxima- 
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tion  is  generally  a  good  one  for  surfaces  with  unknown  characteristics.  Vallado  also  implements  a 
macro-model  approximation  technique  for  the  shape  and  orientation  of  the  satellite  [32].  Vallado’s 
model  will  be  incorporated  into  this  thesis. 

2.3.3  Modified  Equinoctial  Orbital  Elements.  For  most  astrodynamics  problems,  one  of 
the  first  decisions  to  make  is  which  coordinate  frame  to  use  in  defining  the  spacecraft  position  and 
dynamics.  For  the  purposes  of  this  thesis,  geosynchronous  (and  near-geosynchronous)  satellites 
are  being  considered  exclusively.  Due  to  the  special  nature  of  these  orbits,  including  eccentricities 
and  inclinations  of  near  zero  degrees  [4,32],  defining  an  appropriate  coordinate  system  becomes 
important  in  avoiding  singularities.  Perhaps  the  most  obvious  set  of  elements  is  the  equinoctial 
(EQU)  orbital  elements,  as  they  are  free  from  singularities  for  zero  eccentricities  and  both  zero  and 
ninety  degree  inclinations. 

The  EQU  set  are  functions  of  the  mean  anomaly,  M  and  of  the  classical  orbital  elements 
eccentricity,  e,  inclination,  i,  argument  of  perigee,  u>  and  right  ascension  of  the  ascending  node,  Cl. 
Broucke  and  Cefola  applied  the  EQU  elements,  based  on  the  initial  definitions  by  Arsenault  et  ah, 
to  the  two-body  problem  for  general  and  special  perturbations  [6]. 

Somewhat  more  recently,  however,  Walker  et  al.  redefined  the  set  as  the  modified  equinoc¬ 
tial  elements  (MEE)  by  substituting  the  semiparameter  for  the  senrinrajor  axis  and  introducing  a 
different  “fast  variable”,  the  phase  angle.  The  authors  chose  to  replace  the  mean  longitude  of  the 
EQU  set  (Cefola’s  “fast  variable”)  with  the  true  longitude,  thus  setting  the  true  longtitude  as  the 
element  fixing  position  within  the  orbit.  By  replacing  the  semimajor  axis,  the  authors  effectively 
obtained  a  set  of  elements  that  is  applicable  to  most  orbits  and  have  non-singular  equations  of 
motion  [34].  Noting  that  two  of  the  elements  are  singular  for  an  inclination  of  180°  (an  orbit  which 
to  this  date  has  no  satellites  in  it),  the  MEE  set  of  orbital  elements  are  the  most  appropriate  for 
the  application  of  this  thesis  [4] . 
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Several  authors  have  used  the  MEE  set  of  elements  for  projects  and  research  in  the  recent 
years.  Betts  used  the  MEE  set  to  define  both  his  vehicle  dynamics  and  trajectory  modeling  for 
a  low-thrust  optimization  problem  in  2000  [4],  in  a  1994  paper  describing  the  direct  transcription 
method  of  numerical  analysis  for  optimal  interplanetary  low-thrust  orbit  transfers  [3]  and  in  a  2003 
paper  regarding  optimal  low-thrust  trajectories  to  the  moon  [5].  Betts  based  much  of  his  work 
using  the  MEE  element  set  on  theory  developed  by  Dr.  Jean  Kechichian,  who  authored  several 
papers  that  used  EQU  coordinates  [17,18].  The  low-thrust  application  is  of  interest  concerning  this 
thesis,  as  the  acceleration  from  SRP  qualifies,  certainly,  as  a  low-tlrrust  acceleration.  Trawny  et 
al.  used  the  MEE  set  for  a  controlled  lunar  impact  mission,  which  also  included  a  trade  study  of 
low-thrust  pulsed  plasma  thrusters  versus  a  solid  rocket  booster  [30]. 

2.3.4  Solar  Sails.  Much  of  the  recent  emphasis  on  the  exploration  of  SRP  effects  has 
been  in  solar  sail  research.  Refer  to  Figures  2.2  and  2.3  for  representative  illustrations  of  concept 
solar  sails.  Scientists  are  currently  developing  methods  for  harnessing  the  SRP  forces  as  a  source 
of  propulsion  for  long-term  interplanetary  space  travel.  Solar  sails  are  particularly  attractive  for 
this  application  because,  unlike  conventional  propulsion  methods,  solar  sails  require  no  onboard 
fuel  source.  Even  though  the  thrust  produced  from  radiation  pressure  is  small,  it  will  continue  as 
long  as  the  sail  is  viable  and  the  satellite  can  maintain  a  visual  path  to  the  Sun.  The  principles 
that  have  been  developed  for  use  by  solar  sails,  however,  still  apply  to  satellites  not  equipped  with 
sails.  Solar  panels,  temperature-control  panels  and  sun  shades  can  also,  in  effect,  be  used  like  solar 
sails.  This  thesis  will  study  the  feasibility  of  implementing  various  architectures  on  satellites  for 
better  utilization  of  the  perturbing  acceleration  caused  by  SRP  to  raise  a  satellite’s  orbits. 

Victoria  Coverstone  has  spent  a  number  of  years  researching  solar  sails,  and  has  published 
many  papers  regarding  interplanetary  travel  and  transfer  orbits  using  solar  sail  technology.  One 
paper  in  particular  develops  a  technique  for  escaping  the  Earth  using  a  solar  sail.  She  begins 
with  a  spacecraft  in  a  geosynchronous  transfer  orbit  and  uses  a  control  algorithm  to  maximize 
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Figure  2.2:  NASA  Artist’s  Rendering  of  a  Conceptual  Solar  Sail  Deploying  [23]. 
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Figure  2.3: 


Artist’s  Rendering  of  a  Deployed  Circular  Solar  Sail  [23]. 
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the  component  of  sail  force  along  the  velocity  vector,  thus  maximizing  the  instantaneous  rate  of 
increase  in  the  total  orbital  energy.  This  thesis  will  make  use  of  her  work  with  solar  sail  control 
algorithms  in  the  MEE  coordinate  frame  to  optimize  the  force  on  the  satellite,  thereby  raising  the 
orbit  at  a  faster  rate  [11]. 
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III.  Methodology 


The  theory  and  equations  presented  in  this  chapter  were  used  by  the  author  to  build  the  orbital 
propagation  model  in  Matlab®  (reference  Appendix  C  for  the  full  set  of  source  code).  The 
algorithm  to  calculate  the  eccentric  anomaly  was  obtained  from  an  outside  source  (see  Section  C.l). 

3.1  Finding  Classical  Orbital  Elements  from  the  Two-Line  Element  Set 

Data  from  the  most  recent  two-line  element  set  (TLE),  along  with  the  necessary  constants 
for  orbital  propagation,  are  the  required  inputs  into  the  Matlab®  code  in  Appendix  C.  A  TLE 
is  the  standard  format  for  communicating  the  information  pertaining  to  the  orbit  of  a  specific 
Earth-orbiting  object.  A  snapshot  (at  a  given  time,  or  epoch)  of  the  satellite  can  be  taken  and 
sent  to  another  organization  for  purposes  of  tracking  the  object.  More  information  on  the  TLE, 
particularly  for  interpreting  the  TLE,  can  be  found  in  Appendix  A. 

Many  orbital  problems  begin  with  the  calculation  of  the  six  classical  orbital  elements  (COE), 
which  in  most  texts  are  given  by  the  semimajor  axis,  a,  the  eccentricity,  e,  the  inclination,  i, 
the  right  ascension  of  the  ascending  node,  f l  or  RAAN,  the  argument  of  perigee,  to ,  and  the  true 
anomaly,  v.  Figure  3.1  illustrates  the  COE  arrangement  for  a  sample  orbit. 

From  the  TLE,  several  of  the  COE,  along  with  numerous  other  parameters,  are  extracted:  ia 
(degrees),  fl0  (degrees),  e0  (dimensionless),  lo0  (degrees),  mean  anomaly,  Ma  (degrees),  and  mean 
motion,  n0  (revs/day)  [28].  The  subscript  o  indicates  the  value  at  epoch  (reference  time),  which  is 
defined  as  the  date  and  time  given  in  the  TLE. 
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Figure  3.1:  COE  Geometry.  The  COE  a  and  e  are  difficult  to  depict  in  this  figure,  and  thus  are 

omitted. 
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After  converting  degrees  to  radians  and  n  to  units  of  radians  per  second,  the  remaining  COE 


are  calculated  through  the  formulas 


v  —  2  arctan 


(3.1) 

(3.2) 


where  /z©  is  the  gravitational  parameter  of  the  Earth  and  E  is  the  eccentric  anomaly  [14].  Note 
from  this  point  forward  the  subscript  o  will  be  dropped,  as  these  formulas  are  applicable  to  any 
time  t,  in  the  orbit.  E  may  be  calculated  through  a  Newton-Raphson  iteration  (see  Section  C.l  for 
details)  [26]. 

Several  other  orbital  parameters  are  of  interest,  and  may  be  calculated  using  the  equations 


a  (1  +  e) 

(3.3) 

a  (1  —  e) 

(3.4) 

ra  ~  R © 

(3.5) 

rp  -  -R© 

(3.6) 

„  [a? 

2n\  /  — 

(3.7) 

V  m© 

^M© 

2a 

(3.8) 

a  (1  —  ecos  E) 

(3.9) 

r~  R® 

(3.10) 

where  ra  is  the  radius  of  apoapsis,  rp  is  the  radius  of  periapsis,  ha  is  the  height  of  apoapsis,  hp  is 
the  height  of  periapsis,  T  is  the  orbital  period,  e  is  the  total  specific  mechanical  energy,  r  is  the 
instantaneous  radius  of  the  spacecraft  and  h  is  the  instantaneous  height  of  spacecraft  [14]. 
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These  parameters  are  used  to  calculate  the  modified  equinoctial  elements,  and  hence  the 


equations  of  motion,  in  Section  3.2  below. 

3.2  Modified  Equinoctial  Elements 

The  MEE  were  derived  from  a  set  of  equinoctial  orbital  elements  first  described  in  Broucke 
and  Cefola  [6].  The  equinoctial  set  was  derived  to  avoid  the  singularities  that  accompany  the 
COE,  namely,  Q  becoming  indeterminate  as  the  inclination  goes  toward  zero  and  ui  becoming 
indeterminate  as  the  eccentricity  goes  toward  zero.  The  equinoctial  elements,  however  can  result 
in  singular  equations  of  motion  with  some  perturbation  techniques. 

The  MEE  set  was  developed  by  Walker  et  al.  to  eliminate  the  singularities  associated  with 
those  perturbation  techniques.  Walker  et  al.  employed  a  different  “fast  variable” ,  the  true  longitude 
L ,  and  the  semilatus  rectum  p  as  opposed  to  the  mean  longitude  A0  and  the  senrinrajor  axis.  In 
doing  so,  the  authors  derived  a  set  of  elements  that,  when  used  in  any  perturbation  technique, 
enables  them  to  be  applicable  to  all  orbits  and  have  non-singular  equations  of  motion.  These 
characteristics  lend  the  MEE  set  toward  use  for  geosynchronous  orbit  problems.  The  set  of  MEE 
equations  are  defined  in  terms  of  the  COE  as  [34] 


p  =  a  (l  —  e2) 

(3.11) 

/  =  e  cos  (w  +  0) 

(3.12) 

g  =  esin  (w  +  0) 

(3.13) 

h  =  tan  (i/2)  cosf 2 

(3.14) 

k  =  tan  (i/2)  sin  0 

(3.15) 

L  —  0  -(-  oj  T  v 

(3.16) 
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where 


p  =  semiparameter  or  semilatus  rectum 
L  =  true  longitude 

The  usual  terminology  also  defines  vo  =  O  +  u>  as  the  longitude  of  periapsis  and  6  =  u>  +  v  as  the 
argument  of  latitude  [3] . 

Equations  3.11-3.16  can  be  solved  for  the  classical  orbital  elements  in  terms  of  the  modified 
equinoctial  elements  to  facilitate  converting  between  the  two  element  sets,  as  below 


p 

°  1  -P-92 

(3.17) 

e  =  V  f 2  +  g2 

(3.18) 

i  =  2  arctan  (^\/h2  +  k2^j 

(3.19) 

=  arctan  ^2\/h2  +  fc2, 1  —  h2  —  k2^j 

(3.20) 

u>  =  arctan  (g/  f)  —  arctan  (k/h) 

(3.21) 

=  arctan  ( gh  —  fk,  fh  +  gk ) 

(3.22) 

kl  =  arctan  (A;,  h) 

(3.23) 

v  =  L  —  (fl  +  w) 

(3.24) 

=  L  —  arctan  (g/f) 

(3.25) 

U  =  OJ  +  V 

(3.26) 

=  arctan  {h  sin  L  —  k  cos  L,  h  cos  L  +  k  sin  L) 

(3.27) 

where  u  is  the  argument  of  latitude  [3]  used  in  lieu  of  v  for  circular  orbits.  In  the  equations  above, 
expressions  of  the  form  arctan(a,  b)  indicate  a  four  quadrant  inverse  tangent  calculation,  such  as 
atan2  in  Matlab®. 
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In  order  to  implement  the  MEE  set  in  the  two-body  equations  of  motion,  Battin  derived 
the  Cartesian  position  and  velocity  vectors  in  MEE  coordinates.  The  MEE  set  is  related  to  the 
earth-centered  inertial  (ECI),  or  Cartesian,  radius  and  velocity  vectors  through  the  expressions  [2] 


(cos  L  +  a2  cos  L  +  2 hk  sin  L ) 
p-  (sin  L  —  a2  sin  L  +  2 hk  cos  L ) 
( h  sin  L  —  k  cos  L) 


(3.28) 


1  /  £t 

P 


1  IE 


sin  L  +  a2  sin  L  —  2 hk  cos  L  +  g  —  2  fhk  +  a2g ) 
—  cos  L  +  a2  cos  L  +  2 hk  sin  L  —  f  +  2 ghk  +  a2f ) 


2 

77 


j  ( h  cos  L  +  k  sin  L  +  fh  +  gk) 
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where 
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II 
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(3.30) 

s2  =  1  +  h2  +  k2 

(3.31) 

w  =  1  +  f  cos  L  +  g  sin  L 

(3.32) 

r  -  — 

(3.33) 

These  equations  yield  the  Cartesian  state  vector  (r,v)  =  {rx,ry,rz,vx,vy,vz)  when  supplied 
with  the  equinoctial  coordinates  ( p ,  /,  g ,  h ,  k,  L).  Betts  noted  that  the  inverse  transformation  (from 
Cartesian  to  MEE  form)  can  be  calculated,  though  for  this  case  the  true  longitude  L  can  only 
be  found  within  a  multiple  of  27r.  It  must  be  resolved  into  an  angle  between  0  and  2n  through 
knowledge  of  the  reference  epoch  time  [3] . 

Equations  3.11-3.16  are  direct  equations;  they  can  be  individually  differentiated  with  respect 
to  time  to  obtain  the  differential  equations  required  for  solving  orbital  problems  (in  terms  of  the 
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COE  and  their  derivatives).  Walker  et  al.  used  Lagrange’s  planetary  equations 


da  —2a2  dR 

dt  g  dr 

de  —  a(l  —  e)2  dR  1  / 1  —  e2  dR 
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dfl  1  dR 
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(3.34) 

(3.35) 

(3.36) 

(3.37) 

(3.38) 


and,  expressing  the  classical  elements  in  terms  of  the  modified  equinoctial  elements  and  performing 
the  appropriate  transformation  of  R  (a,  e,  i,  w,ft,t  —  r)  into  R  ( p ,  /,  g ,  h,  k.  L),  derived  the  Gaussian 
equations  of  motion 
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k  =  a - - - sin  L 

y  g  2w 

ij  =  '/pp(—\  + —a  (hsinL  —  kcosL)  Ar 


\P 


w  y  fi 


(3.39) 

(3.40) 

(3.41) 

(3.42) 

(3.43) 

(3.44) 


where  Ar,  Af  and  A„  are  the  non- two-body  perturbation  accelerations  in  the  radial,  tangential 
and  normal  directions,  respectively.  The  convention  for  the  radial  direction  is  along  the  geocentric 
radius  vector  of  the  spacecraft  (measured  positive  in  the  direction  away  from  the  geocenter);  the 
tangential  direction  is  perpendicular  to  the  radius  vector  (measured  positive  in  the  direction  of 
orbital  motion);  the  normal  direction  is  the  positive  direction  along  the  angular  momentum  vector 
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of  the  satellite’s  orbit  [3,34],  as  defined  in  Equations  3.49-3.51  below.  Betts  mentions  in  his  paper 
that  Equation  3.41  reflects  the  derivation  performed  by  Battin  [2]  and  corrects  an  error  presented 
in  the  original  paper  by  Walker  et  al.  [3,34]. 


3.3  Equations  of  Motion 


The  modified  equinoctial  dynamics  equations  can  be  rearranged  into  matrix  form  [4],  ex¬ 


pressed  as 

y  =  AA  +  b 


(3.45) 


where 


A  = 


£  sin  L 
n 


j-  cos  L 


2p  p 

w  y  p 


p_L 

p  w 


[(to  +  1)  cosi  + /]  - 


^  [(«'  +  !)  sin  L  +  g\ 


-  —  (h  sin  L  —  k  cos  L) 

p  W  v  J 


-  —  (h  sin  L  —  k  cos  L) 

LL  W  '  ' 


p_  s2  cos  L 
p  2  w 


p  2  w 


^  (hsinL  —  kcosL) 


(3.46) 


and 


b  = 


0  0  0  0  0 


(3.47) 


The  total  disturbing  acceleration  vector  A  is  given  by 


A  —  A  rir  T  A  tit  -f-  A  nin 


(3.48) 
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which  is  expressed  in  a  rotating  radial  coordinate  frame  whose  principle  axes  are  defined  by 


%f  — 


r 


(3.49) 


in 


f  X  V 

||f  x  v\\ 


it  —  in  x  ir 


(r  x  v)  x  r 
||f  x  u||||f|| 


(3.50) 

(3.51) 


These  are  the  unit  vectors  in  the  radial,  normal  and  tangential  directions,  respectively.  The  3x3  ro¬ 
tation  matrix,  R,  from  the  radial-tangential-normal  (RTN)  coordinate  frame  to  the  ECI  coordinate 
frame  is  defined  in  terms  of  the  unit  position  vectors  in  Equations  3.49-3.51  as 


R  = 


Xf 


it 


in 


(3.52) 


For  unperturbed  two-body  motion,  A  is  zero.  The  first  five  equations  of  motion  then  become 
p  =  f  —  g  =  h  =  k  =  0,  indicating  that  all  of  the  modified  equinoctial  elements  except  for  the  fast 
variable,  L,  are  constant.  It  is  also  important  to  note  that  any  perturbing  forces  can  be  included 
in  A,  given  that  they  are  rotated  (if  necessary)  into  the  RTN  coordinate  frame.  For  a  satellite  in 
geosynchronous  orbit,  the  three  most  influential  perturbations  are  those  caused  by  Earth  oblateness, 
gravitational  disturbances  from  secondary  bodies  and  SRP  [32] . 

3.3.1  Gravitational  Disturbing  Acceleration.  Earth  oblateness,  which  is  sometimes  re¬ 
ferred  to  as  simply  the  J2  perturbation  when  higher  order  zonal  harmonics  (as  well  as  all  sectorial 
and  tesseral  harmonics)  are  neglected  (see  Figure  3.2),  induces  periodic  variations  in  all  orbital 
elements  (a,  e,  i,  f 2,  u>,  M).  Secular  perturbations,  however,  are  only  induced  by  even  zonal  har¬ 
monics,  and  only  affect  Cl,  to  and  M  [32],  For  the  purposes  of  this  thesis,  only  the  J2  perturbations 
will  be  modeled,  as  these  have  the  most  effect  on  a  satellite  within  a  geosynchronous  orbit  about 
the  Earth.  The  models  often  used  for  determining  the  perturbing  acceleration  due  to  oblate  body 
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Figure  3.2:  Disturbing  Forces  on  the  Gravity  Field  of  a  Central  Body.  The  shaded  regions  in 

subfigures  (a),  (b)  and  (c)  indicate  areas  of  increased  gravity. 

(a)  Zonal  Harmonics.  Zonal  harmonics  often  coincide  with  the  more  massive  latitudinal  bands 
of  the  Earth  (e.g.,  the  equator).  They  account  for  most  of  the  non-spherical  Earth  gravitational 
anomalies. 

(b)  Sectorial  Harmonics.  The  sectorial  harmonics,  so  called  due  to  the  resemblance  to  sections  of 
an  orange,  correlate  to  longitudinal  bands  of  increased  mass  (e.g.,  North-South  running  mountain 
ranges  like  the  Rockies  and  Appalachians). 

(c)  Tesseral  Harmonics.  These  harmonics  represent  specific  regions  on  the  Earth  which  depart  from 
the  perfect  sphere  model,  e.g.,  Australia  and  Antartica. 


effects  are  normally  defined  in  a  local  horizontal  reference  frame.  The  non-spherical  gravitational 
acceleration  vector  is  given  in  this  frame  as 


Sg 


SgNiN  -  8grir 


(3.53) 


where  defines  the  unit  vector  in  the  local  North  direction 


with 


£n 


t 


0  0  1 


(3.54) 


(3.55) 


and  ir  defines  the  unit  radial  direction  component  as  given  in  Equation  3.49. 

Ignoring  the  tesseral  and  sectorial  harmonics,  the  general  forms  for  the  oblate  Earth  pertur¬ 
bations  to  the  gravitational  acceleration  model  (in  the  North  and  radial  directions,  respectively) 
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are 


SgN 


Sgr 


p'kJk 

(3.56) 

Pk  J k 

(3.57) 

where  /z®  is  the  Earth’s  gravitational  constant,  r  is  the  geocentric  distance  to  the  satellite,  P® 
is  the  equatorial  radius  of  the  Earth,  (f>gc  is  the  geocentric  latitude,  Jk  is  the  fcth  zonal  harmonic 
coefficient  and  Pk  (sin (f>gc)  is  the  fcth  order  Legendre  polynomial  with  corresponding  derivative  P(.. 
For  a  model  that  only  includes  zonal  harmonics,  the  east  component  of  the  zonal  gravity  effects  is 
exactly  zero  and  can  be  set  as  such  to  generate  a  fairly  accurate  model  [4,5]. 

For  J2  effects  only,  the  perturbations  to  the  gravitational  acceleration  are 


SgN 

5gr 


cos 

?,2 


p^j2 


where  Vallado  [32]  defines 


(3.58) 

(3.59) 


J2  =  0.0010826269  (3.60) 

P2  =  i(3sin20gc-l)  (3.61) 


The  derivative  of  the  Legendre  polynomial  is  (note  that  there  is  also  a  cos  4>gc  term  that  has  already 
been  absorbed  in  Equation  3.58) 

P2  =  3sin</>ffc  (3.62) 
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and  the  geocentric  latitude  is  easily  computed  using  the  ECI  position  vector  of  the  satellite 


4>gc 


(3.63) 


where  the  /,  J  and  K  components  of  the  radius  vector  are  given  as  77,  rj  and  vk  respectively. 
The  gravitational  perturbation  in  the  rotating  radial  reference  frame,  then,  is  achieved  through  a 
rotation  by  matrix  R 

A g  =  R  TSg  (3.64) 


3.3.2  Solar  Radiation  Pressure  Perturbing  Acceleration.  Like  atmospheric  drag,  SRP  is 
a  nonconservative  non-gravitational  disturbance.  The  SRP  perturbing  force  increases  in  relative 
magnitude,  however,  as  the  effects  of  atmospheric  drag  decrease.  In  GEO,  for  example,  the  perturb¬ 
ing  force  due  to  atmospheric  drag  can  be  considered  to  be  negligible,  whereas  SRP  is  the  dominant 
perturbation  at  that  altitude. 

SRP  perturbing  acceleration  is  difficult  to  model  for  many  reasons.  Most  models  implement 
a  constant  solar  flux,  though  solar  flux  in  fact  varies  both  throughout  the  year  and  during  periods 
of  solar  activity.  The  variations  associated  with  periods  of  increased  solar  activity  cannot  be 
accurately  forecasted  due  to  their  sporadic  tendencies,  but  the  variations  exhibited  in  solar  flux  are 
computationally  insignificant  for  the  purposes  of  this  thesis  (on  the  order  of  7%)  due  to  the  larger 
ambiguities  in  surface  characteristics. 

For  an  accurate  model,  the  cross-sectional  areas  of  the  satellite  must  also  be  computed  with 
reasonable  precision,  and  the  Earth’s  shadowing  effects  on  the  satellite  must  be  considered.  Other 
difficulties  arise  in  estimating  the  reflective  properties  of  the  satellite,  as  the  composition  of  the 
satellite  surfaces  will  directly  affect  how  much  force  from  SRP  can  be  harnessed  to  produce  incre¬ 
mental  acceleration.  The  development  of  an  accurate  model  also  involves  the  ability  to  calculate 
the  precise  position  of  the  Sun  with  respect  to  the  satellite  and  the  Earth  at  any  given  time  during 
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the  orbit,  as  well  as  the  correct  attitude  of  the  satellite  (and  orbit)  and  the  exact  value  of  the 
SRP  itself.  Using  this  determination,  the  satellite  can  be  reoriented  with  respect  to  the  Sun  at  the 
proper  positions  to  maximize  (and  in  turn  minimize)  the  components  of  the  force  vector  in  the  unit 
normal  direction. 


To  account  for  the  yearly  variance  of  the  solar  flux,  Wertz  developed  a  time-varying  formula 
to  replace  the  often-used  solar  flux  constant  SF  =  1353  W/m2.  His  equation  for  the  solar  flux, 
below,  uses  a  value  measured  from  the  Earth  aphelion  point,  which  can  be  obtained  from  tables 
in  the  Astronomical  Almanac.  An  abbreviated  form  of  the  tables  is  presented  in  Appendix  B  [31]. 
D aphelion  is  defined  as  2it  times  the  number  of  days  elapsed  from  Earth  aphelion,  as  a  fraction  of 
the  whole  year. 


SF  = 


_ 1358 _ 

1.004  +  0.0334  COS  D aphelion 


w_ 

m2 


(3.65) 


The  change  in  momentum,  or  the  force  of  the  solar  pressure  per  unit  area  psr,  is  defined  in 
terms  of  the  solar  flux  as 


SF 

PSR  = - 

c 


(3.66) 


where  c  is  the  speed  of  light. 

Most  models  for  SRP  perturbations  do  not  include  the  complex  surface  geometry  of  a  satellite. 
Satellite  surfaces  degrade  with  time  while  in  orbit,  and  although  reflective  properties  may  have  been 
known  with  relative  accuracy  at  launch,  over  many  years  the  coefficients  will  fluctuate  drastically 
and  models  must  be  updated  to  reflect  the  true  perturbing  forces  experienced  by  the  satellite. 
To  incorporate  a  more  accurate  picture  of  how  the  satellite  behaves  under  perturbing  forces  from 
SRP,  an  algorithm  must  include  both  the  absorption  and  reflection  of  a  representative  solar  ray. 
Figure  3.3  illustrates  the  geometry  involved  when  calculating  the  forces  on  the  satellite  due  to  SRP. 

The  two  primary  types  of  reflection,  specular  and  diffuse,  are  illustrated  in  Figure  3.4.  The 
solar  radiation  incident  on  the  satellite  surface  area,  as  in  Figure  3.3,  shows  the  incident  angle  <t>inc , 
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Figure  3.3:  Incident  Solar  Radiation  Geometry.  Note  that  h  and  in  are  oriented  upwards  from 

the  plane  of  the  paper;  s  and  rsat@  are  collinear  with  the  incident  solar  ray.  The  unit  vectors  t 
and  n  are  the  unit  vectors  in  the  tangential  and  normal  directions  (with  respect  to  the  flat  plate), 
respectively. 
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defined  as  the  angle  between  the  normal  unit  vector  h  and  the  unit  vector  defining  the  direction 
from  the  satellite  to  the  Sun,  s.  Note  that  the  reflected  angle,  (j)ref ,  and  the  incident  angle  are 
equal  for  purely  specular  reflection  [35] . 

Vallado  derives  a  relationship  for  the  perturbing  acceleration  by  assuming  a  Lambertian 
diffusion  (see  Section  2.3.2)  to  model  the  diffuse  and  specular  radiation  forces  as  they  affect  the 
satellite.  He  uses  a  macro-model  to  sum  each  satellite  surface  area,  as  follows, 

-  _  VsrA 0i  cos  <j>  inci  \  0  | CRdi  ^  l^i/i  \~'l  fo  a>-7  \ 

dSR  =  -)J - -12  — — I-  CRSi  cos  (pina  n+(l  —  CRsi)s>  (3.67) 

z '  m  l  o  J  ) 

i= 1 

where  A&i  is  the  average  effective  cross-sectional  area  of  surface  i,  m  is  the  mass  of  the  entire 
satellite,  and  and  CRSi  are  the  coefficients  of  diffuse  and  specular  reflectivity,  respectively  [32] . 
Note  that  the  coefficient  of  absorption  is  not  explicitly  needed,  although  it  has  been  incorporated 
since  1  —  Crs  =  Cro,  +  Cud-  The  surface  normal  unit  vector  h  and  the  solar  incidence  unit 
vector  s  are  required  to  specify  the  orientation  of  the  satellite.  A  macro-model  estimates  the 
shape  and  orientation  of  the  satellite  to  find  a  more  accurate  acceleration;  the  summation  adds 
the  contributions  from  the  n  flat  plate  surfaces  of  the  satellite  model,  where  each  surface  is  the  ith 
element. 

For  the  case  of  a  spin-stabilized  satellite  with  one  rotating  platform  and  one  despun  platform, 
such  as  the  DSCS  II  satellite,  Equation  3.67  simplifies  to  (with  n  —  2) 

asR  =  ~asRd  ~  asRa  (3.68) 

where  d  and  s  indicate  the  despun  and  spin-stabilized  platforms,  respectively,  and 

& SRd  —  PSRAexpd  COS  (j)inCd  \ ^  0  T  CRsi  COS  (f)incd  Rd  T  (1  C-Rsd)  ^(3  (3.69) 

to  t  L  3  J  J 

RSRs  ~  ~^PSRAexps  COS  (j)inca  —  b  CRs,  COS  4>incs  Rs  T  (1  CRss )  (3.70) 
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(a) 


Diffuse  surface 


(b) 

Figure  3.4:  (a)  Specular  Reflection.  A  mirror-like  reflection,  specular  reflection  is  often  identified 

by  the  congruent  angles  the  incoming  and  reflected  waves  make  with  the  surface  normal  vector, 
(b)  Diffuse  Reflection.  With  this  type  of  reflection,  the  wave  bounces  off  of  the  reflecting  surface 
at  seemingly  random  angles,  due  to  the  surface  irregularities. 
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The  average  effective  cross-sectional  areas,  Aexpa  and  Aexpdl  are  approximated  by  flat  plates  (refer 
to  Figure  3.5).  For  cylinders,  this  is  a  valid  approximation,  as  the  area  as  seen  from  the  Sun  will 
indeed  be  represented  as  a  flat  plate  area.  For  the  parabolic  antennas  on  the  DSCS  II  satellite’s 
despun  platform,  flat  plate  approximations  are  used  due  to  the  ambiguity  of  the  measurements 
available. 


Figure  3.5:  DSCS  II  F-13  Approximate  Area  Illustration.  Although  the  despun  platform  ap¬ 

proximate  area  is  shown  with  some  depth  in  the  figure,  for  the  purposes  of  this  thesis  it  is  assumed 
to  have  zero  depth.  The  spin-stabilized  platform  is  faced  with  solar  panels. 

Since  the  focus  of  the  thesis  is  to  raise  the  orbit  of  a  geosynchronous  satellite  the  most 
effectively  (read:  most  quickly)  with  only  SRP,  it  is  also  interesting  to  consider  larger  (at  least 
conceptually)  flat  plate  areas  for  the  DSCS  II  parabolic  antennae  approximation.  The  motivation 
for  this  thesis  is  satellite  disposal,  and  it  is  unlikely  that  very  large  solar  sail-like  devices  would 
be  used  for  this  purpose.  Therefore,  only  smaller  areas  are  considered;  the  original  flat  plate 
approximated  area  (area  factor  equal  to  1)  will  be  multiplied  by  two  different  area  factors:  five  and 
ten.  The  simulation  will  be  repeated  for  each  of  these  areas,  and  then  the  change  in  semimajor  axis 
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will  be  plotted  versus  area  factor  (1,  5  and  10),  as  well  as  the  time  to  disposal  versus  area  factor. 
These  plots  should  give  an  indication  of  how  the  satellite  responds  to  the  SRP  perturbation  with 
changing  flat  plate  areas. 


The  coefficients  of  absorption  and  diffuse  and  specular  reflectivities  are  also  considered  un¬ 
known  values.  The  satellite  used  in  the  model  has  been  in  orbit  since  1979,  and  the  since  the 
satellite  surfaces  have  degraded  over  time,  the  current  coefficients  cannot  be  accurately  estimated. 
For  the  purposes  of  this  thesis,  the  coefficients  (note  that  CRa  +  cr d  +  crs  =  1)  have  each  been  set 
to  one-third  unless  noted.  Several  other  test  cases  will  also  be  simulated:  all  absorption  (cRa  =  1), 
all  specular  reflection  (crs),  all  diffuse  reflection  ( CRd )  and  a  blend  of  specular  and  diffuse  reflection 
(cRd  =  0.5  and  crs  =  0.5).  This  will  highlight  the  contribution  of  these  types  of  reflection  to  the 
orbit  raising  solution. 


Specific  parameters  regarding  the  DSCS  II  satellite,  including  the  TLE  set  used  for  orbit 
propagation,  can  be  found  in  Appendix  A  [23]. 


The  unit  vector  defining  the  direction  from  the  satellite  to  the  Sun  at  any  given  time,  s,  can 
be  computed  through 


^  satQ 

ll^satoll 


(3.71) 


where  rsato  is  the  vector  from  the  satellite  to  the  Sun.  From  Figure  3.6,  it  is  evident  that  rsato  is 
simply  the  sum  of  two  vectors 

rSate  =  Abo  -  f  (3.72) 


where  r®®  is  the  Earth-Sun  vector  and  f  is  the  ECI  position  vector  of  the  satellite. 

The  Earth-Sun  vector,  at  any  given  time  t,  is  calculated  using  an  algorithm  from  Vallado’s 
text.  Given  the  current  Julian  date,  the  number  of  Julian  centuries  is  calculated  as 


JD  —  2451545.0 
TDB  ~  36525 


(3.73) 
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Satellite 


Figure  3.6:  Earth-Satellite-Sun  Vector  Diagram.  The  bodies  and  distances  in  the  figure  are  not 
to  scale;  they  are  positioned  for  purposes  of  clarity.  The  plane  of  the  ecliptic  is  assumed  congruent 
with  the  plane  of  the  paper. 

where  Ttdb  references  the  number  of  Julian  centuries  based  on  barycentric  dynamical  time  (TDB). 

Using  the  number  of  Julian  centuries,  the  mean  longitude  of  the  Sun  A mq  ,  solar  mean  anomaly 
M0,  ecliptic  longitude  A ec;,  distance  between  the  Earth  and  the  Sun  rg0  and  the  obliquity  of  the 
ecliptic  e  are  calculated  from  the  equations 


M0  =  357.5277233°  +  35999.05034  TTDb  (3.74) 

Am0  =  280.460°  +  36000.770  TTDb  (3.75) 

A  ed  =  A  m0  +  1.914666471  sin  (M0)  +  0.019994643  sin  (2M0)  (3.76) 

re0  =  1.000140612-  0.016708617  cos  (M0)  -  0.000139589  cos  (2Me)  (3.77) 

e  =  23.439291  —  0.0130042  Ttdb  (3.78) 


The  ecliptic  latitude  of  the  Sun  is  never  greater  than  0.000333°,  and  is  set  to  zero  for  this  problem 
without  much  loss  of  accuracy.  The  Earth-Sun  vector,  in  Astronomical  Units  (AU),  is 


COS  (A eel) 


Ubo  —  re© 


cos(e)  sin(Aeci) 
sin(e)  sin(Aeci) 


(3.79) 
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To  convert  from  AU  to  kilometers,  multiply  Equation  3.79  by  149598000  [32]. 

The  parameters  for  the  satellite’s  rotating  platform  to  use  in  Equation  3.70  are  easily  com¬ 
puted.  The  DSCS  II’s  spinning  platform  is  of  cylindrical  shape;  the  reflected  or  incident  angle 
of  the  cylinder  is  always  <j>s  =  0  since  the  face  of  the  cylinder  will  always  be  perpendicular  to  s. 
Subsequently,  the  unit  normal  vector  h  will  always  lie  along  s. 

The  computation  of  the  unit  normal  vector  for  the  despun  platform  is  more  complicated. 
Because  the  parabolic  antennae  are  not  cylindrical,  as  the  satellite  changes  its  attitude,  the  effec¬ 
tive  cross-sectional  area  presented  to  the  Sun  also  changes.  The  control  algorithm  illustrated  in 
Section  3.4  outlines  several  methods  for  calculating  n  for  the  parabolic  antennae. 

3.3.3  Solar  Gravitational  Perturbations.  For  a  satellite  orbiting  the  Earth,  the  magnitudes 
of  the  perturbations  caused  by  the  effects  of  gravitational  attraction  from  a  third  body  (e.g.,  Sun 
or  Moon)  become  greater  as  the  altitude  of  the  orbit  increases.  The  third-body  effects,  as  they  are 
termed,  become  noticeable  as  the  effects  of  atmospheric  drag  subside;  they  are,  in  fact,  the  largest 
perturbations  affecting  satellites  in  geosynchronous  orbit.  The  primary  third-body  gravitational 
perturbing  forces  acting  on  a  satellite  in  geosynchronous  orbit  are  those  from  the  Sun  and  Moon, 
although  the  Moon’s  perturbations  will  not  be  considered  due  to  the  complexity  of  the  calculations 
and  the  associated  processing  time.  Since  the  Moon  does  not  induce  secular  changes  in  senrinrajor 
axis  of  the  orbit,  this  exclusion  does  not  diminish  the  results. 

In  a  geosynchronous  orbit,  the  disturbances  caused  by  the  Sun  can  have  a  significant  long-term 
perturbing  effect  on  the  orbit.  Vallado  lists  the  main  effect  as  a  regression  of  the  satellite’s  orbital 
plane  about  an  axis  normal  to  the  Sun’s  orbital  plane.  The  only  secular  perturbations  caused  by 
the  solar  gravitational  perturbation  appear  in  the  node  and  perigee.  The  expected  long-periodic 
variations  in  e,  i,  O  and  lo  are  entirely  due  to  the  regression  of  the  satellite’s  perigee  and  the  motion 
of  the  Earth  about  the  Sun.  When  modelled  perfectly,  the  Sun  causes  no  secular,  long-periodic 
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or  rn-monthly  variations  in  semimajor  axis,  the  parameter  of  most  interest  in  this  thesis  [32].  To 
generate  a  more  complete  model,  however,  these  third-body  perturbations  will  be  included. 

The  equations  used  to  calculate  the  perturbing  acceleration  from  the  Sun  are  nearly  identical 
to  the  previous  section.  The  equation  of  motion  for  a  three-body  system  (in  this  section,  the 
Earth-Sun-satellite  system)  is  given  by 


-  +^0 

'  ®sa£ 


(  T satQ 
\  ^ satO 


where  the  perturbing  acceleration  due  to  the  Sun  is  isolated  as 


(3.80) 
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and  where  //©  is  the  gravitational  parameter  of  the  Sun,  fsatQ  is  the  position  vector  from  the 
satellite  to  the  Sun  and  r©©  is  the  position  vector  from  the  Earth  to  the  Sun  (the  calculation  of 
which  is  presented  in  Equations  3.74-3.79). 

Although  at  this  point  the  equation  could  be  numerically  integrated,  closer  study  yields  a 
potential  problem.  The  distances  from  the  satellite  to  the  Sun  and  the  Earth  to  the  Sun  are 
very  similar  and  thus  the  value  in  the  parenthesis  will  be  very  small.  It  is  possible  that  this 
could  induce  errors  during  a  numerical  integration  of  the  equation,  unless  the  capability  exists  to 
store  large  numbers  of  significant  figures  (which  would  significantly  lengthen  the  time  required  for 
integration) .  A  series  expansion  of  r~atQ  yields  a  more  numerically  stable  solution  of 

7^  M©  ( —  o-.  r^sat  ■  r®d  /r®sat-f®0\  \  cox 

d(lQ  =  --3—  r®sat  -  3re© - 3 - 15re©  I  - 5 -  1  (3.82) 
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where  an  expansion  of  three  terms  for  the  solar  perturbation  case  is  sufficiently  accurate  [32] . 
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The  perturbing  acceleration  in  the  rotating  radial  reference  frame  is  [5] 


A90=RT 


5qe 


(3.83) 


3-4  Control  Algorithm 

If  left  to  drift  uncontrolled,  a  satellite  in  geosynchronous  orbit  will  experience  only  periodic 
motion  in  the  semimajor  axis  due  to  SRP  perturbations  [32],  This  is  due  to  the  Earth’s  rotation 
around  the  Sun;  for  one-half  of  the  orbit,  SRP  acts  to  add  a  Av  to  the  perigee  of  the  satellite’s  orbit 
and  subtract  a  Av  from  the  apogee,  thus  causing  the  orbit  to  gradually  become  elliptic.  For  the 
remaining  six  months  of  Earth’s  orbit,  SRP  causes  the  satellite  to  experience  the  opposite  reaction: 
a  Av  is  subtracted  from  the  perigee  and  added  to  the  apogee,  circularizing  the  satellite’s  orbit 
again  (see  Figure  3.7  for  further  clarification).  The  purpose  of  this  thesis,  however,  is  to  harness 
the  SRP  perturbing  acceleration  and  raise  the  orbit  of  a  satellite  into  an  appropriate  disposal  orbit 
(the  guidelines  of  which  can  be  found  in  Appendix  D).  To  raise  the  orbit  of  a  satellite  the  most 
effectively  with  respect  to  time,  a  host  of  different  methods  could  be  implemented.  Two  methods 
will  be  addressed,  the  first  of  which  is  a  controller  that  would  only  require  attitude  changes  twice 
daily.  The  second  is  a  more  complex  algorithm  that  continuously  changes  the  satellite’s  attitude. 
Both  methods  require  only  electrical  power  to  be  implemented,  which  for  the  DSCS  II  is  provided 
by  the  available  solar  arrays  mounted  on  the  spin-stabilized  platform. 


3-4-1  Simple  Controller.  Perhaps  the  simplest  control  algorithm  is  one  based  on  the 
trigonometries  of  the  orbital  problem.  Using  the  diagram  in  Figure  3.8,  it  is  evident  that  the  cosine 
of  the  angle  9  can  be  computed  using  the  equation 


cos  9  = 


T' satQ  ’  VSat 
||  T sat©  ||  ||  ^sat  || 


(3.84) 
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Figure  3.7:  The  Behavior  of  a  Geosynchronous  Satellite  Over  One  Year  Under  the  Long-Periodic 
Changes  in  Semimajor  Axis  Due  to  SRP. 

If  cos  8  is  negative,  the  satellite  is  moving  away  from  the  Sun;  if  it  is  positive,  the  satellite 
is  moving  toward  the  Sun.  The  control  algorithm  computes  cos  9 ,  and  if  negative,  positions  the 
parabolic  antennae  perpendicular  to  re®,  thus  maximizing  the  SRP  perturbation  addition  to  the 
acceleration  and  raising  the  orbit.  This  sets  the  effective  area  to  its  maximum  value  (for  DSCS 
II,  approximately  3.3528  m2).  If  positive,  the  controller  positions  the  antenna  platform  parallel  to 
r®®  to  minimize  the  SRP  perturbation  contribution.  This,  in  effect,  “turns  off”  the  SRP  for  half 
of  the  satellite’s  orbit,  permitting  the  orbit  raising  to  dominate  instead  of  being  cancelled  out  by 
the  orbit  lowering.  In  the  code,  the  thickness  of  the  flat  plate  area  is  approximated  to  be  0  m;  this 
sets  the  antenna  area  to  0  m2 . 

The  main  advantage  to  a  simplified  control  algorithm  is  less  operator  intervention.  For  an 
older  satellite  already  in  orbit,  such  as  DSCS  II,  there  is  no  realistic  opportunity  to  install  a  complex 
controller  or  to  reprogram  the  onboard  computer.  Thus,  feedback  from  the  operators  suggests  that 
a  satellite  operator  would  manually  need  to  change  the  satellite’s  attitude  for  every  re-orientation 
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Figure  3.8:  Simple  Controller  Diagram.  Note  that  this  diagram  is  not  to  scale.  Cosine  is  negative 
above  the  horizontal  axis,  which  is  defined  by  t^q,  and  is  positive  below  the  horizontal  axis. 

required  (certainly  tedious).  This  twice-daily  algorithm  is  one  that  would  result  in  the  least  amount 
of  operator  hours  spent  raising  the  orbit. 

The  downside  of  the  simple  controller  is  directly  related  to  its  advantageous  feature  of  less 
operator  oversight.  If  changing  the  attitude  of  the  satellite  twice  per  orbit  gleans  a  moderate 
change  in  semimajor  axis,  it  follows  that  the  more  reorientations  performed  over  one  orbit,  the  more 
resulting  change  in  semimajor  axis.  For  this  reason,  a  more  complex  controller  will  be  studied,  one 
that  constantly  reorients  the  satellite  to  maximize  the  change  in  semimajor  axis. 

3-4-2  Complex  Controller.  For  newer  satellites,  or  satellites  yet  to  be  launched  or  devel¬ 
oped,  there  are  many  other  possible  control  algorithms  that  would  be  able  to  raise  the  orbit  of  a 
satellite  more  quickly  than  the  simple  algorithm  in  Section  3.4.1. 

Coverstone  and  Prussing  [11]  derived  a  sail  force  algorithm  that  forces  the  unit  normal  vector, 
n,  to  maximize  the  rate  of  orbital  energy,  E.  Their  controller  continuously  orients  the  sail  (in  three 
dimensions)  to  maximize  the  component  of  the  sail  force  along  the  satellite’s  velocity  vector,  v.  This, 
in  turn,  maximizes  the  instantaneous  rate  of  increase  of  the  total  orbital  energy.  The  authors’  full 
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derivation  can  be  found  in  their  journal  article;  however,  it  is  worth  discussing  the  highlights  of  the 
algorithm  and  the  changes  implemented  for  this  thesis. 

Coverstone  and  Pressing  made  several  simplifying  assumptions  in  their  algorithm  that  con¬ 
tradict  the  desired  model  for  this  thesis.  They  derive  a  simplified  perturbing  acceleration  model 
in  which  they  assume  a  constant  solar  flux  and  only  specular  reflection.  They  introduce  an  ac¬ 
celeration  vector  that  has  no  component  in  the  satellite-Sun  direction,  indicating  that  they  have 
neglected  the  forces  due  to  diffuse  reflection  and  absorption.  The  simplified  acceleration  is  given  as 


^  2  WA  T  , 

b  = - (n  a) 

me 

=  So  (nTd)  n 


(3.85) 


where  W  is  the  constant  solar  intensity  (flux)  at  1  AU,  A  is  the  sail  area,  c  is  the  speed  of  light, 
h  is  the  unit  normal  vector  to  the  sail  and  a  is  the  unit  vector  in  the  direction  of  the  Sun.  Using 
Equation  3.85,  the  authors  derive  the  equation  for  the  time  rate  of  change  of  the  orbital  energy  as 


E  = 


(3.86) 


where  <Jq  is  the  gravitational  perturbation  of  the  Sun.  The  authors  note,  at  this  point  in  their 
derivation,  that  their  solution  is  not  precisely  the  minimum  time  solution,  although  for  low  acceler¬ 
ation  rates  (as  in  this  thesis,  with  a  very  small  “sail” )  their  solution  approaches  the  minimum  time 
solution.  They  concede  that  for  many  applications,  a  shorter-duration  solution  may  be  achieved, 
however,  since  no  propellant  is  being  used,  having  the  real  minimum-time  trajectory  is  probably 
not  a  critical  issue. 

To  maximize  the  orbital  energy,  the  authors  orient  the  sail  such  that  the  expression  STv 
is  maximized  in  Equation  3.86  above.  In  other  words,  this  will  maximize  the  component  of  the 
sail  acceleration  along  the  instantaneous  velocity  vector  at  each  point  in  the  trajectory.  In  the 
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Figure  3.9:  a(3Z  Axis  Definitions.  Note  that  Coverstone  et  al.  assume  the  heliocentric  Z-axis  is 

collinear  with  the  a(3Z  Z-axis. 


derivation,  the  components  of  the  unit  normal  vector  h  that  result  in  the  maximization  of  E,  in  an 
a/3Z  reference  frame  (see  Figure  3.9),  are 


r  _  -M 

(3.87) 

np  =  £,na 

(3.88) 
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where 
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and  va  and  vp  are  rotated  into  the  a/3z  reference  frame  through  a  rotation 


(3.90) 


Va 

cos  *F  sin  ’F  0 

vx 

Vp 

= 

—  sin  4/  cos  'F  0 

Vy 

Vz 

0  0  1 

Vz 

(3.91) 


where  \F  is  the  celestial  longitude  of  the  Sun  [11], 
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The  algorithm  used  in  this  thesis  differs  from  this  point  forward.  Coverstone  and  Prussing, 
in  deriving  their  algorithm,  neglected  the  effect  of  the  inclination  of  the  Earth’s  axis  with  respect 
to  the  ecliptic  plane  (approximately  23.5°).  For  the  purposes  of  this  thesis,  an  additional  rotation 
is  performed  between  the  ECI  and  heliocentric  (XYZ)  reference  frames.  Using  the  ecliptic  latitude 
from  Equation  3.76,  the  rotation  matrix  from  the  ECI  frame  to  the  XYZ  frame  is 


R 


XYZECI 


I  0  0 

0  cos  23.5  sin  23.5 
0  —sin  23.5  cos  23.5 


Then,  v  in  the  heliocentric  frame  is 


(3.92) 
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and  in  the  a/3Z  frame, 
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(3.94) 


Once  the  components  of  the  unit  normal  vector  in  the  a/3Z  frame  are  computed  using  Equa¬ 
tions  3.87-3.89,  the  desired  h  is  then  rotated  back  to  the  ECI  frame  using  the  transposes  of  the 
matrices  in  Equations  3.92  and  3.94.  Note  that  to  implement  the  authors’  equations  with  Vallado’s 
algorithm  for  the  perturbing  SRP  acceleration,  the  negative  of  Coverstone  and  Prussing’s  unit 
normal  vector  must  be  used.  This  puts  the  measurement  of  4>inc  in  the  same  frame  as  Vallado, 
between  the  unit  normal  n  and  the  unit  vector  in  the  satellite-Sun  direction  s  versus  between  the 
acceleration  vector  in  Equation  3.85  and  the  negative  of  the  satellite-Sun  vector. 
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The  appropriate  <f>inc  for  the  despun  platform  is  calculated  using  the  following  relationship 


<fiinc  =  arccos  (n  •  s)  (3.95) 

and  the  calculations  for  the  SRP  perturbing  accelerations  resume  at  Equation  3.68. 

3.5  Earth  Shadow  Effects 

Over  long  periods  of  time,  it  has  been  determined  that  Earth’s  shadow  plays  a  secular  role  in 
the  perturbation  of  a  satellite’s  orbit  [7,10].  As  stated  in  Chapter  II,  the  inclusion  of  such  a  shadow 
function  in  this  thesis  is  particularly  necessary,  primarily  due  to  the  desire  to  raise  the  senrinrajor 
axis  of  an  orbit.  Bryant  declared,  in  1961,  that  the  addition  of  an  accurate  shadow  function  to  an 
orbital  propagation  algorithm  would  in  fact  result  in  a  (positive)  secular  change  in  semimajor  axis. 
Without  the  model  of  Earth’s  shadow,  the  SRP  perturbing  acceleration  would  only  direct  periodic 
changes  in  the  senrinrajor  axis  [7]. 

Two  particular  categories  of  shadow  functions  exist:  a  simplified  cylindrical  model,  as  in 
Ferraz  Mello’s  model  [12],  and  a  conical  model,  such  as  that  presented  by  Kabelac.  While  the  conical 
method  is  more  accurate,  it  is  often  difficult  to  implement  [16].  Kabelac’s  model,  in  particular, 
uses  a  more  realistic  conical  Earth  shadow  with  two  main  regions,  the  penumbra  and  umbra  (see 
Figure  3.10  [32]).  The  umbra  is  the  region  that  is  in  total  eclipse  by  the  Earth  (or  other  primary 
body  for  applications  with  any  planetary  body),  and  the  penumbra  encompasses  the  region  only 
partially  eclipsed  by  the  Earth  or  other  primary  body.  The  model  also  accounts  for  the  mutual 
position  of  the  Sun,  Earth  and  satellite,  the  shape  of  conical  surface,  light  diffusion,  atmospheric 
absorption  and  the  effect  of  ozone  [16]. 

Ferraz  Mello’s  model,  however,  approximates  the  Earth’s  shadow  as  a  cylinder,  and  hence  has 
only  one  region  (see  Figure  3.11  [32]).  This  approximation  steins  from  an  assumption  that  the  Sun 
is  at  an  infinite  distance  from  the  Earth.  His  algorithm,  while  greatly  simplified  from  Kabelac’s, 
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Satellite 


Figure  3.10:  Eclipse  Geometry  with  Penumbra  and  Umbra  Regions.  The  penumbra  is  the  lighter 
gray  region  and  the  umbra  is  the  dark  gray  region.  The  figure  is  looking  down  onto  the  ecliptic 
plane;  it  is  not  to  scale. 


Figure  3.11:  Cylindrical  Eclipse  Geometry.  The  Sun  is  approximated  as  a  point  source,  hence 

the  shadow  cast  by  the  Earth  is  cylindrical  in  shape  instead  of  conical.  The  figure  is  not  to  scale, 
and  the  satellite  orbit  is  not  intended  to  approximate  the  orbit  of  DSCS  II. 


enables  the  user  to  include  a  shadow  model  without  significant  loss  of  accuracy.  Kabelac  states,  in 
his  paper,  that  his  model  differs  only  on  the  order  of  0.1%  to  1%  in  the  orbital  elements  [12, 16].  A 
version  of  Ferraz  Mello’s  model  is  used  in  this  thesis. 

The  portion  of  the  orbit  when  a  satellite  is  fully  inside  the  cylindrical  shadow  region  of  the 
Earth  is  of  most  interest.  When  an  orbiting  object  is  in  the  Earth’s  shadow,  the  SRP  perturbing 
acceleration  has  no  effect  on  the  orbit  of  the  satellite;  the  eclipsed  satellite  is  not  exposed  to  SRP 
due  to  obstruction  from  the  Earth.  Vallado  derived  an  algorithm  to  determine  if  two  objects  have 
line  of  sight,  and  this  can  be  directly  applied  to  the  shadow  problem.  The  algorithm  can  be  used  to 
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Figure  3.12:  Geometry  of  Tmin.  Line  of  Sight  exists  between  the  two  positions  a  and  b  if  the 

sum  of  the  two  angles  is  larger  than  the  calculated  value  for  6  (the  angle  between  the  two  vectors). 
The  parametric  approach,  in  which  the  location  for  the  minimum  distance  to  the  Earth  is  found, 
can  also  be  used;  each  vector  can  be  represented  as  a  function  of  r  [32], 


determine  if  the  satellite  has  a  direct  line  of  sight  with  the  Sun,  and  hence  whether  it  is  in  Earth’s 
shadow.  Vallado  assumes  in  his  derivation  that  the  light  from  the  Sun  acts  as  a  point  source.  He 
derives  a  minimized  parametric  value,  Tmj„,  that  minimizes  the  distance  to  the  central  body  (see 
Figure  3.12  for  an  illustration  of  the  geometry  involved  in  the  definition  of  rmjn) 


7~min 


Ini2  -fi-f2 
\n\2  +  \r2\2  -  2ri  ■  f2 


(3.96) 


Next,  he  derives  a  parametric  representation  of  a  line  between  the  two  position  vectors  fj  and  r2, 
and  substitutes  the  minimum  parametric  value  Tmin  to  get 


|c(Trmn)|  — 


2  (1  -  Tmin)  |ri  |2  +  (fl  •  f2)  Tn 


Rl 


(3.97) 


The  radius  of  the  Earth  squared,  172. ,  in  km,  is  present  in  the  denominator  to  ensure  (if  r)  and  f2 
are  in  units  of  km)  that  the  value  of  the  parametric  function  is  in  units  of  Earth  radians  [32] .  Im¬ 
plementing  the  algorithm  in  the  orbit  propagation  code  requires  an  if  statement  within  Matlab®. 
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A  summarized  version  is  presented  below;  the  complete  algorithm  can  be  found  in  Appendix  C. 


If  T~min  ^  0-0  Of  Tmin  ^1-0 

the  satellite  is  illuminated 
elseif  |c(rmi„)|2  >  1.0 

the  satellite  is  illuminated 

else 


the  satellite  is  not  illuminated 


(3.98) 


If  the  satellite  is  illuminated,  the  SRP  perturbing  acceleration  algorithm  is  applied;  conversely,  if 
the  satellite  is  not  illuminated,  the  SRP  perturbing  acceleration  vector  is  set  to  zero. 
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IV.  Results 


Each  set  of  results  in  this  chapter  includes  one  of  three  controllers:  no  controller,  the  simple  control 
algorithm  or  the  complex  controller.  All  results  include  a  cylindrical  Earth  shadow  model,  and  the 
perturbations  caused  by  third-body  gravitational  acceleration,  oblate  Earth  (J2)  and  SRP.  They 
incorporate  Vallado’s  most  sophisticated  SRP  perturbation  algorithm  with  varying  solar  flux  and 
the  inclusion  of  absorption  and  reflection  coefficients  (Equation  3.67).  Three  individual  areas  were 
considered  for  the  despun  platform:  the  original  area  (reference  Section  3.3.2  for  a  description  of 
how  this  value  was  estimated),  the  original  area  multiplied  by  a  factor  of  five  and  the  original  area 
multiplied  by  a  factor  of  ten.  Unless  explicitly  noted,  the  results  (a  tabular  summary  of  which  can 
be  found  in  Section  4.7)  presented  use  the  original  area.  It  is  also  important  to  note  that  in  all  the 
plots  throughout  this  chapter,  “time”  on  the  x-axis  indicates  the  time  (in  seconds)  elapsed  since 
epoch,  as  defined  in  Section  3.1. 

Jf.,1  Baseline  Case 

The  baseline  case,  one  with  no  active  controller  to  reorient  the  satellite  to  maximize  SRP 
effects,  is  an  interesting  one  for  purposes  of  comparison.  Essentially,  the  satellite  is  drifting  around 
the  Earth  freely  in  a  nominal  state,  that  is,  without  manual  reorientation.  Thus,  one  would  expect 
the  satellite  to  maintain  a  circular  orbit  without  secular  changes  in  senrinrajor  axis,  much  like  the 
two-body  problem.  Additionally,  the  perturbations  should  generate  only  periodic  variations  in  e 
and  i,  and  both  secular  and  periodic  changes  in  u>,  fl  and  v. 

The  following  plots  are  the  results  for  the  DSCS  II  F-13  satellite  with  no  controller,  i.e.  the 
apparent  area  of  the  despun  platform  is  set  to  zero,  and  the  projected  area  of  the  spin-stabilized 
platform  is  set  to  its  maximum  of  5.0168  m2 . 

Validation  results  for  the  baseline  cases  (of  1,  100  and  1000  orbits),  propagated  with  Satellite 
Tool  Kit  (STK),  are  also  presented  in  this  section.  This  effort  was  requested  by  SMC  Det  12  to 
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determine  if  the  Matlab®  code  developed  by  the  author  yielded  reasonable  results  as  compared 
to  a  “known”  solution  from  trusted  source  (e.g.  STK).  The  other  cases  in  the  thesis  were  unable 
to  be  generated  in  STK  for  purposes  of  validation,  due  primarily  to  the  use  of  controllers. 

The  area  factor  referenced  in  Section  3.3.2  only  affects  the  flat-plate  area  (despun  platform 
apparent  area);  the  results  for  varying  areas  will  yield  the  same  results  as  the  original  despun 
platform  area,  as  this  is  set  to  zero  for  all  the  baseline  cases.  Note  also  that  any  changes  in  values 
for  the  absorption  and  specular  and  diffuse  reflection  coefficients  will  not  affect  the  solution,  as 
these  are  also  multiplied  by  the  despun  platform  area  (zero  for  the  baseline  case)  in  Equation  3.67. 
To  avoid  duplication,  therefore,  only  one  set  of  plots  will  be  generated  for  each  of  the  three  baseline 
cases,  and  no  further  baseline  plots  will  be  included  in  the  sections  to  come. 

To  verify  that  the  thesis  results  for  the  baseline  cases  are  indeed  reasonable,  the  satellite’s 
motion  was  propagated  in  STK  for  each  case.  For  purposes  of  comparison,  the  HPOP  propagator 
was  used  to  simulate  the  orbit  with  the  same  perturbations  as  considered  in  the  thesis.  For  the  Earth 
gravity,  the  model  JGM2.grv  was  chosen,  with  a  maximum  degree  of  2  and  a  maximum  order  of  0. 
Solar  radiation  pressure  was  implemented  with  a  coefficient  of  reflection  of  Cr  =  1  (corresponding 
to  all  specular  reflection),  an  area-to-mass  ratio  of  0.00821604  m2 /kg,  and  a  cylindrical  shadow 
model.  The  gravitational  perturbations  due  to  the  Sun  were  also  included.  Atmospheric  drag, 
tides  and  all  other  gravitational  bodies  were  excluded. 

Figures  4. 1-4. 5  highlight  the  characteristics  of  the  DSCS  II’s  first  orbit  from  epoch.  Fig¬ 
ures  4.1  and  4.2  depict  the  behavior  of  the  ECI  position  and  velocity  vectors  as  the  satellite  orbits 
once  around  the  Earth.  Figure  4.1  illustrates  the  periodic  motion  of  the  x,  y  and  3  components  of 
the  ECI  position  and  velocity  vectors.  As  would  be  expected  for  a  circular  orbit,  the  2  components 
exhibit  very  little  change. 

Although  the  two-dimensional  orbit  trace  in  Figure  4.2  appears  circular  at  first  glance,  it  is  in 
fact  slightly  elliptical  (eG  =  0.00056).  The  graph  is  a  plot  of  of  the  x  and  y  components  of  the  ECI 
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Position  vector  versus  time  over  1  orbital  periods 


Velocity  vector  versus  time  over  1  orbital  periods 


Figure  4.1:  ECI  Position  and  Velocity  Vectors  Versus  Time  over  One  Orbit. 


2-D  orbit  trace  of  satellite  over  1  orbital  periods 


Figure  4.2:  The  Two-dimensional  Orbit  Trace  over  One  Orbit.  The  starting  position  of  the 

satellite  (at  epoch)  is  indicated  by  a  filled  circle  labeled  by  “Starting  value”  on  the  right  side  of  the 
graph. 
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position  vector,  and  as  such  is  “looking  down”  on  the  equatorial  plane  of  the  Earth  (represented 
as  the  filled  circle  in  the  center  of  the  figure). 

Figures  4.3  and  4.4,  depict  the  behavior  of  i,  a  and  e  over  the  corresponding  time  period. 
The  inclination  of  the  initial  orbit  oscillates  around  the  epoch  value  of  14.4995°.  For  satellites  in 
a  geosynchronous  orbit,  i  drifts  to  about  15°  over  the  course  of  about  27  years  [32],  so  one  would 
expect  that  an  inclination  in  the  range  of  0°  to  15°  is  a  realistic  value.  The  periodic  motion  in  i 
is  common;  perturbations  due  to  zonal  harmonics,  drag,  third-bodies  and  SRP  all  cause  short  or 
long  periodic  tenancies  in  inclination  [32] . 

In  this  thesis,  all  secular  changes  in  the  satellite’s  semimajor  axis  are  a  direct  result  of  the 
controller  usage.  For  low  Earth  orbit,  conversely,  the  drag  perturbations  induce  secular  changes 
in  a.  The  changes  in  semimajor  axis  (Aa)  throughout  the  chapter,  excluding  the  Aa  indicated 
in  Figures  4.3,  4.6  and  4.8,  which  are  simply  the  mean  values  over  the  integration  time,  are  cal¬ 
culated  using  the  polyfit  function  in  Matlab®.  It  is  important  to  note  that  the  polyfit  function 
becomes  more  accurate  as  the  number  of  periodic  cycles  increases,  and  hence  the  values  given  by 
the  function  for  one  orbit  are  not  necessarily  the  most  accurate  representation  of  the  change  in 
semimajor  axis.  The  apparent  change  exhibited  in  Figure  4.3  is  a  result  of  the  short  period  (high 
frequency)  oscillations  (from  the  J2,  third-body  and  SRP  perturbations  [32])  affecting  the  mean 
value  calculation.  Similarly,  the  apparent  secular  change  in  eccentricity  is  solely  a  factor  of  the 
short-  and  long-periodic  variations  caused  by  those  perturbations.  No  additional  plots  of  e,  i,  u>,  fl 
or  v  will  be  included  in  beyond  this  section,  as  the  semimajor  axis  is  the  parameter  of  interest. 

Figure  4.3  illustrates  one  orbit  propagated  with  STK  overlaid  on  one  orbit  from  the  Matlab® 
code  results.  The  STK  results  are  very  similar  to  the  thesis  results  as  would  be  expected,  however, 
the  inherent  differences  found  between  the  method  STK  uses  to  propagate  an  orbit  and  the  methods 
presented  in  this  thesis  will  yield  slightly  different  answers.  The  mean  values  of  the  satellite’s 
semimajor  axis  are  labeled  on  the  first  subplot  for  each  method.  These  are  a  serviceable  indication 
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of  the  comparative  methods,  but  the  more  interesting  statistics,  the  difference  between  the  vectors 


at  each  time  step  (deviation)  are  given  in  subplot  2. 


x  i  o4  Semimajor  axis  versus  time 


Deviation  of  STK  and  Matlab  methods 


Figure  4.3:  STK  Similar  Case  Validation  Results  for  Baseline  Case  Behavior  of  a  for  1  Orbit. 


Figure  4.5  depicts  the  behavior  of  u>,  Q  and  v  for  one  orbit.  The  STK  results  are  plotted 
with  the  Matlab®  results  for  purposes  of  comparison.  Again,  the  different  methods  yield  very 
similar  results  for  the  osculating  elements.  The  periodic  motion  of  ui  is  expected;  all  perturbations 
from  zonal  harmonics,  drag,  third-bodies,  and  SRP  induce  periodic  motion.  Any  secular  motion, 
although  small,  should  be  expected  as  well;  even  zonal  harmonics,  third-body  perturbations  and 
SRP  perturbations  affect  the  motion  of  the  satellite  secularly  [32].  fl  and  i/,  while  uninteresting 
for  such  a  short  period,  do  reflect  the  predictable  behavior  of  these  COE  for  the  satellite,  and  thus 
indicate  an  accurate  algorithm.  The  precession  of  v  over  the  course  of  the  orbit  is  evident  in  the 
figure.  The  right  ascension,  f2,  exhibits  a  downward  secular  trend,  indicating  that  the  ascending 
node  (the  point  on  the  equatorial  plane  where  the  satellite  crosses  South  to  North)  is  precessing 
toward  the  Vernal  Equinox  (in  a  westerly  direction). 
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Figure  4.4: 
Orbit. 
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Figure  4.5:  STK  Similar  Case  Validation  Results  for  Baseline  Case  Behavior  of  w,  f 1  and  v  for 

1  Orbit. 
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Figures  4.6  and  4.7  illustrate  the  results  for  100  orbits  from  STK  and  Matlab®.  Again, 
note  that  the  STK  results  are  very  close  to  the  Matlab®  results.  The  increase  in  deviation 
over  time  is  most  likely  indicative  of  small  time  step  differences  in  the  respective  integrators  of 
STK  and  Matlab®  (HPOP  and  ode45),  not  of  a  growing  difference  between  the  semimajor  axis 
vectors  over  the  course  of  many  orbits.  Since,  especially  for  the  100  orbit  case,  the  deviation 
appears  to  have  nearly  equal  discrete  steps,  another  explanation  for  the  rising  difference  stems  from 
Matlab® Matlab®  using  double  precision  and  STK  using  a  single  precision-type  integrator.  If 
looked  at  closely,  the  trace  of  a  for  both  methods  lies  nearly  colinearly  over  the  entire  time  set. 


x  -|  o4  Semimajor  axis  versus  time 


Deviation  of  STK  and  Matlab  methods 


Figure  4.6:  STK  Similar  Case  Validation  Results  for  Baseline  Case  Behavior  of  a  for  100  Orbits. 

Figures  4.8  and  4.9  illustrate  the  results  for  1000  orbits  from  STK  and  Matlab®.  As  in  the 
previous  plots  of  1  and  100  orbits,  the  semimajor  axis  exhibits  nearly  zero  change,  as  would  be 
expected  under  no  controller  influence.  The  Matlab®  results  are  again  similar,  in  both  period 
and  magnitude,  to  the  STK  results.  Notice  that  the  mean  deviation  has  decreased  from  the  100 
orbit  case,  and  the  1  sigma  standard  deviation  value  has  as  well. 
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Inclination  versus  time 


time  (sec) 


Figure  4.7:  STK  Similar  Case  Validation  Results  for  Baseline  Case  Behavior  of  i  and  e  for  100 

Orbits. 
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Figure  4.8:  STK  Similar  Case  Validation  Results  for  Baseline  Case  Behavior  of  a  for  1000  Orbits. 
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Inclination  versus  time 


Figure  4.9:  STK  Similar  Case  Validation  Results  for  Baseline  Case  Behavior  of  i  and  e  for  1000 

Orbits. 

In  conclusion,  note  that  the  results  generated  exhibited  only  periodic  changes  in  semimajor 
axis,  inclination  and  eccentricity,  as  expected.  The  periodic  and  secular  variations  were  present  in 
the  plots  of  w,  f 1  and  v,  again  as  expected.  Just  as  notable,  the  STK  results  seems  to  validate 
the  Matlab®  code;  the  small  differences  between  the  methods  are  most  likely  due  to  only  the 
numerical  integration  differences  between  HPOP  and  ode45.  Most  importantly,  both  methods 
induced  no  secular  changes  in  semimajor  axis  over  any  number  of  orbits. 


4-2  Nominal  Coefficient  Case 

This  section  addresses  the  cases  in  which  the  satellite’s  coefficients  of  absorption,  specular 
reflection  and  diffuse  reflection  are  each  set  to  equal  one-third  (hereafter  referred  to  as  the  nominal 
coefficient  case).  The  nominal  case  was  initially  chosen  to  represent  the  current  properties  of 
the  DSCS  II  satellite  since  they  are  not  well  known.  Recall  that  the  surfaces  of  the  satellite 
have  degraded  significantly  since  launch  into  orbit,  and  the  current  coefficients  of  absorption  and 
reflection  cannot  be  directly  measured.  Without  good  data  on  the  coefficients,  equal  weighting 


58 


is  a  reasonable  starting  point  for  analysis.  Later  in  the  chapter,  other  possible  values  of  the 
reflection  and  absorption  coefficients  will  be  investigated  and  used  in  the  orbital  propagation  of 
the  satellite,  including  the  cases  where  each  coefficient  is  made  to  be  dominant.  In  particular,  four 
additional  cases  will  be  examined  to  observe  the  changing  effects  on  the  satellite’s  perturbed  orbit: 
all  absorption  (cRa  =  1),  all  specular  reflection  (crs  =  1),  all  diffuse  reflection  ( end  =  1)  and  a  mix 
of  diffuse  and  specular  reflection  ( crs  =  0.5  and  cud  =  0.5). 

For  this  case,  and  the  cases  to  follow,  one  would  expect  secular  changes  in  semimajor  axis  due 
to  the  use  of  the  controllers.  Since  the  area-to-mass  ratio  factors  linearly  into  Equation  ??  for  the 
perturbing  acceleration  due  to  SRP,  one  would  also  expect  that  the  higher  the  ratio,  the  greater 
the  secular  change  in  semimajor  axis  (with  a  linear  relationship).  In  general,  one  would  expect 
the  complex  controller  to  yield  greater  secular  changes  in  semimajor  axis,  but  this  could  prove  to 
be  dependent  on  the  designated  coefficients  of  absorption  and  reflection  since  the  controller  was 
designed  for  purely  specular  reflection. 

4-2.1  Simple  Controller  Case.  This  section  covers  the  test  case  results  for  the  simple 
controller  described  in  Section  3.4.1.  Again,  the  nominal  coefficients  of  reflection  and  absorption 
were  used.  Two  subcases  were  considered:  one  using  the  original  approximated  despun  platform 
area  and  the  second  using  varying  areas  to  illustrate  the  effects  of  larger  exposed  areas  to  the  Sun. 
Note  that  each  of  the  remaining  plots  for  1000  orbits  has  been  orbit  averaged,  that  is,  the  data 
sets  have  been  reduced  via  the  Matlab®  polyfit  command  to  a  linear  best  fit  curve  so  as  to  give 
a  clearer  representation  of  the  results. 

4-2. 1.1  Original  Despun  Platform  Area.  The  original  approximated  area  of  the 
satellite’s  despun  platform  (based  on  the  cross-sectional  areas  of  the  parabolic  antennae)  is  3.3445 
m2 .  The  results  in  this  section  include  that  area  in  the  Matlab®  algorithm. 
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For  ease  of  comparison,  the  behavior  of  the  semimajor  axis  using  the  simple  controller  are 
overlaid  on  the  corresponding  baseline  results  from  Section  4.1.  The  behaviors  of  i,  e,  to,  and 
v  are  for  the  most  part  uninteresting.  In  any  case,  as  long  as  the  results  are  as  expected  per  the 
perturbations  included  in  the  algorithm,  their  values  do  not  affect  the  desired  solution  (ultimately, 
the  raise  in  semimajor  axis).  Except  for  the  baseline  cases  in  the  previous  section,  therefore,  the 
plots  of  i  and  e  will  be  omitted  in  favor  of  plots  solely  of  the  semimajor  axis  over  time. 

Figure  4.10  is  a  comparison  of  the  baseline  and  simple  controller  cases  for  one  orbit  only.  Note 
that  in  each  plot  hereafter,  the  change  in  semimajor  axis  for  the  controller  case  (calculated  with 
polyfit)  is  given  in  the  title  above  the  plot.  The  change  in  semimajor  axis,  A  a,  for  the  baseline  case 
is  0.35462  km ;  for  the  simple  controller  case,  it  is  0.30981  km..  Normally,  one  would  expect  that 
using  the  simple  controller  would  be  more  advantageous  than  using  no  controller.  The  rise  over 
one  orbital  period,  however,  is  difficult  to  calculate  accurately  because  the  change  in  senrinrajor 
axis  is  based  on  the  line  fit  of  the  osculating  element  a,  rather  than  a  computation  of  the  mean 
elements.  Essentially,  the  short-periodic  variations  in  a  are  adversely  skewing  the  calculation  of  Aa 
for  one  orbit.  The  line  fit  method  gives  more  representative  results  when  used  for  the  100  and  1000 
orbit  results.  The  plots  for  100  and  1000  orbits,  as  illustrated  in  Figures  4.11  and  4.12  respectively, 
yield  better  results  than  the  baseline  case.  For  100  orbits,  the  secular  change  in  semimajor  axis  is 
3.186  km,  while  for  1000  orbits  it  is  31.6495  km. 

4- 2. 1.2  Varied  Despun  Platform  Area.  The  results  presented  in  this  section  use  two 
area  factors  to  illustrate  the  changes  in  behavior  of  the  semimajor  axis.  The  first  plot  was  generated 
with  the  area  factor  equal  to  five;  the  second  was  generated  with  the  area  factor  set  equal  to  ten. 

The  1000  orbit  results  with  an  area  factor  of  five  (the  original  despun  platform  area  multiplied 
by  five)  are  shown  in  Figure  4.13.  At  this  point,  the  results  begin  to  vary  significantly  from  the 
results  using  the  nominal  parabolic  antennae  effective  area,  especially  for  the  larger  number  of 
orbits.  Although  the  plot  is  not  shown,  for  one  orbit,  the  change  in  semimajor  axis  using  the 
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Figure  4.10:  Simple  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  1  Orbit 

and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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Figure  4.11:  Simple  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  100  Orbits 
and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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Figure  4.12:  Simple  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  1000 

Orbits  and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

simple  controller  is  0.49196  km.  For  100  orbits,  the  gain  from  using  the  simple  controller  over  no 
controller  is  15.8675  km  and  for  1000  orbits  it  is  158.7075  km.  This  is  a  marked  improvement  over 
the  results  in  Section  4. 2. 1.1. 

The  1000  orbit  results  with  an  area  factor  of  ten  (the  original  despun  platform  area  multiplied 
by  ten)  are  shown  in  Figure  4.14.  Once  again,  the  results  vary  significantly  from  the  nominal 
parabolic  antennae  effective  area,  especially  for  the  larger  number  of  orbits.  For  one  orbit,  the 
gain  from  using  the  simple  controller  over  no  controller  is  0.3650  km.  For  100  orbits,  the  change 
in  semimajor  axis  is  31.7532  km  and  for  1000  orbits  it  is  318.2681  km.  Based  on  these  results,  it 
appears  that  the  rate  of  change  of  semimajor  axis  is  directly  proportional  to  the  variable  area,  as 
would  be  expected  from  the  equation  for  the  force  due  to  SRP  (Equation  3.67). 


4-2.2  Complex  Controller  Case.  This  section  covers  the  test  case  results  for  the  complex 
controller  described  in  Section  3.4.2.  Again,  the  nominal  coefficients  of  reflection  and  absorption 
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Figure  4.13:  Simple  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  1000 

Orbits  and  an  Area  Factor  of  5,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


Figure  4.14:  Simple  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  1000 

Orbits  and  an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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were  used.  Two  subcases  were  considered:  one  using  the  original  approximated  despun  platform 
area  and  the  second  using  varying  areas  to  illustrate  the  effects  of  larger  parabolic  antennae. 

4- 2. 2.1  Original  Despun  Platform  Area.  Using  the  satellite’s  original  apparent 
despun  platform  area,  the  semimajor  axis  values  for  the  complex  controller  and  baseline  cases 
for  1000  orbits  were  plotted  against  the  elapsed  time  in  Figure  4.15.  For  one  orbit,  the  raise  in 
semimajor  axis  is  0.31468  km.  For  100  orbits,  the  difference  from  using  the  complex  controller  over 
having  no  controller  is  3.4137  km,  and  for  1000  orbits  it  is  34.2982  km.  The  complex  controller, 
for  the  nominal  coefficient  case,  yields  approximately  an  8%  increase  over  the  simple  controller. 


Figure  4.15:  Complex  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  1000 

Orbits  and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

4-2. 2. 2  Varied  Despun  Platform  Area.  The  results  presented  in  this  section  use  two 
area  factors  (five  and  ten)  to  illustrate  the  changes  in  behavior  of  the  semimajor  axis.  The  plot 
shown  was  generated  with  the  area  factor  set  equal  to  ten. 

Similarly  to  the  results  for  the  simple  controller,  a  larger  effective  area  for  the  despun  platform 
results  in  a  significant  benefit  in  semimajor  axis  increase.  For  an  area  factor  of  five,  over  one  orbit 
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the  semimajor  axis  rose  0.51633  km,  resulting  in  a  benefit  over  having  no  controller  of  0.16171  km). 
Over  100  orbits,  the  benefit  is  17.0713  km,  and  over  1000  orbits  it  is  172.1135  km. 

Figure  4.16  shows  the  1000  orbit  results  when  the  despun  platform  effective  area  is  multiplied 
by  ten.  Over  one  orbit,  the  benefit  over  having  no  controller  is  0.41377  km.  For  100  orbits,  it  is 
34.1493  km,  and  for  1000  orbits  it  is  345.1880  km. 


Figure  4.16:  Complex  Controller  Case  Behavior  of  a  for  Nominal  Reflection  Coefficients,  1000 

Orbits  and  an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

As  expected,  the  semimajor  axis  results  yielded  by  using  the  controllers  were  dominated  by 
the  secular  changes,  although  the  periodic  changes  are  evident  in  the  plots  for  1  and  100  orbits.  For 
at  least  the  nominal  case,  and  it  is  predicted  for  the  remaining  cases,  the  area  factors  do  appear  to 
yield  a  linear  relationship  with  the  changes  in  semimajor  axis. 

4-3  Absorption  Case 

This  section  details  the  results  from  the  cases  in  which  the  diffuse  and  specular  coefficients 
of  reflection  have  been  set  to  equal  zero  and  the  coefficient  of  absorption  has  been  set  to  one 
(hereafter  referred  to  as  the  “all-absorption”  case).  It  is  expected  that  the  performance  of  the 
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satellite  (with  regards  to  raising  the  semimajor  axis)  under  the  SRP  perturbation  will  be  the  least 
in  the  all-absorption  case,  primarily  due  to  the  direction  of  the  resultant  force  vector  under  the 
all-absorption  condition. 

4-3.1  Simple  Controller  Case.  This  section  addresses  the  results  for  the  all- absorption 
case  using  the  simple  controller  described  in  Section  3.4.1.  Two  subcases  were  considered:  one 
using  the  original  approximated  despun  platform  area  and  the  second  using  varying  areas  (of  five 
and  ten  times  the  original  area)  to  illustrate  the  effects  of  a  larger  despun  platform  area  on  the 
semimajor  axis. 


4-3. 1.1  Original  Despun  Platform  Area.  The  results  in  this  section  use  an  apparent 
despun  platform  area  of  3.3445  ?n2.  Similarly  to  Section  4.2. 1.1,  the  semimajor  axis  results  from 
the  simple  controller  test  cases  are  overlaid  on  the  baseline  cases  for  purposes  of  comparison. 

In  Figure  4.17,  note  that  the  two  cases  nearly  overlap.  For  one  orbit,  the  raise  in  semimajor 
axis  for  the  original  area  case  is  0.29399  km,  while  for  the  baseline  case  it  is  0.35462  km..  The  all 
absorption  result  is,  as  expected,  slightly  less  than  the  nominal  case  result  in  Section  4.2.1. 

The  data  for  100  and  1000  (see  Figure  4.18)  orbits  yield  similarly  scaled  results.  For  100 
orbits,  the  change  in  semimajor  axis  is  2.1010  km.,  and  for  1000  orbits  it  is  20.7914  km..  These 
results  are,  again,  slightly  lower  the  nominal  case. 

4-3. 1.2  Varied  Despun  Platform  Area.  The  results  in  this  section  use  two  different 
area  multiplication  factors  to  illustrate  the  changes  in  semimajor  axis  when  using  larger  despun 
platform  areas.  The  first  plot,  Figure  4.19,  uses  an  area  factor  of  five  applied  to  the  despun 
platform’s  original  effective  area.  The  second,  Figure  4.20,  uses  an  area  factor  of  ten. 

The  simple  controller  does  not  glean  much  benefit  over  the  no  controller  case  for  one  orbit. 
The  advantage  is  only  0.05823  km.,  compared  to  0.1373  km  for  the  nominal  case.  The  difference 
for  100  orbits  is  10.4413  km  and  for  1000  orbits  (Figure  4.19)  it  is  104.3604  km.  Using  the  nominal 
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Figure  4.17:  Simple  Controller  Case  Behavior  of  a  for  all  Absorption,  1  Orbit  and  an  Area  Factor 
of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


1  q«  a  versus  time  over  1 000  orbital  periods;  A  a  =20.791 4  km 


Figure  4.18:  Simple  Controller  Case  Behavior  of  a  for  all  Absorption,  1000  Orbits  and  an  Area 

Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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coefficients  of  absorption  and  reflection,  A  a  for  1000  orbits  was  158.7181  km ;  for  the  all  absorption 
case,  it  is  104.3709  km.  As  expected,  the  all  absorption  case  yields  a  lower  semimajor  axis  increase 
than  the  nominal  case. 


Figure  4.19:  Simple  Controller  Case  Behavior  of  a  for  all  Absorption,  1000  Orbits  and  an  Area 

Factor  of  5,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

The  resulting  Aa  for  an  area  multiplication  factor  of  ten  are,  of  course,  higher  than  the 
previous  results  for  the  all  absorption  test  case.  The  change  in  semimajor  axis  for  one  orbit  is 
0.56143  km ,  for  100  orbits  is  20.9048  km  and  for  1000  orbits  (Figure  4.20)  is  209.1176  km.  Over 
1000  orbits,  therefore,  the  semimajor  axis  increase  is  about  100  km.  less  than  in  the  corresponding 
nominal  case,  where  it  was  318.2681  km. 

4-3.2  Complex  Controller  Case.  Perhaps  counterintuitively,  the  implementation  of  the 
complex  controller  yields  less  semimajor  axis  change  than  the  simple  controller  for  the  all  absorption 
case,  as  illustrated  in  Figure  4.22.  The  complex  controller  reorients  the  despun  platform  more  often 
than  the  simple  controller,  thus  the  apparent  area  of  the  despun  platform  will  be  greater  for  a  longer 
period  of  time  during  the  orbit.  This  causes  the  resultant  force  vector  from  SRP  to  dominate  the 


68 


Figure  4.20:  Simple  Controller  Case  Behavior  of  a  for  all  Absorption,  1000  Orbits  and  an  Area 

Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

complex  controller  case  and  the  effects  of  all  absorption,  therefore,  degrade  the  performance  of  the 
satellite  under  SRP  perturbations  more  quickly.  In  short,  the  complex  controller  is  optimized  for 
purely  specular  reflection,  not  for  pure  absorption.  The  percent  difference  in  using  the  complex 
controller  versus  the  simple  controller  is  approximately  —39%. 

4-3.2. 1  Original  Despun  Platform  Area.  This  section  addresses  the  results  of  using 
the  complex  controller  for  the  all  absorption  case  with  the  original  despun  platform  approximated 
area.  The  resulting  changes  in  the  semimajor  axis,  although  an  improvement  on  the  baseline  case, 
indicates  that  this  is  the  least  effective  scenario  for  raising  the  semimajor  axis  of  the  satellite  orbit 
(among  the  cases  studied  for  the  purposes  of  this  thesis). 

Over  one  orbit,  A  a  is  0.28844  km,  which  is  0.06618  km  less  than  the  baseline  case.  For  100 
orbits,  A  a  is  1.28454  km.  For  1000  orbits,  as  illustrated  in  Figure  4.21,  A  a  is  12.6063  km.  This 
is  approximately  sixty  percent  of  the  semimajor  axis  raise  from  using  the  simple  controller  (see 
Section  4. 3. 1.1)  and  about  a  third  of  the  complex  controller  nominal  case  value  in  Section  4. 2. 2.1. 
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Figure  4.21:  Complex  Controller  Case  Behavior  of  a  for  all  Absorption,  1000  Orbits  and  an  Area 
Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

4- 3. 2. 2  Varied  Despun  Platform  Area.  The  results  presented  in  this  section  use  the 
two  aforementioned  area  factors  of  five  and  ten.  Figure  4.22  illustrates  the  1000  orbit  results  using 
an  area  factor  of  ten  applied  to  the  despun  platform. 

Like  the  original  area  results  in  Section  4. 3. 2.1,  the  changes  in  semimajor  axis  in  this  section 
continue  to  be  less  than  all  other  cases  except  the  baseline  case.  In  particular,  for  an  area  factor 
of  five  over  one  orbit,  Aa  is  0.38511  km..  For  100  orbits,  A  a  is  6.4357  bn,  and  for  1000  orbits  it 
is  63.2965  km.  For  1000  orbits,  the  raise  in  semimajor  axis  is  about  40  km  less  than  that  for  the 
simple  controller  under  all  absorption. 

Figure  4.22  illustrates  the  effects  of  having  ten  times  the  original  estimated  despun  platform 
area  for  1000  oribts.  For  one  orbit,  the  resulting  Aa  of  0.50595  km  is  only  marginally  better  than 
for  an  area  factor  of  five.  Aa  is  12.8627  km  for  100  orbits,  and  126.7702  km  for  1000  orbits,  which 
is  significantly  less  than  for  the  complementary  nominal  case. 
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Figure  4.22:  Complex  Controller  Case  Behavior  of  a  for  all  Absorption,  1000  Orbits  and  an  Area 
Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

4-4  Specular  Reflection  Case 

The  specular  reflection  section  addresses  the  test  cases  in  which  the  coefficient  of  specular 
reflection  of  the  satellite’s  despun  platform  is  set  to  one  and  the  remaining  coefficients  are  set  to  zero. 
As  in  the  previous  sections,  several  subcases  were  considered.  Test  cases  were  run  in  Matlab® 
for  both  the  simple  controller  and  the  complex  controller;  for  each  controller  the  original  despun 
platform  approximated  area  and  the  area  factors  of  five  and  ten  times  the  original  area  were  used. 


4-4-1  Simple  Controller  Case.  This  section  presents  the  results  for  the  purely  specular 
reflection  case  using  the  simple  controller  described  in  Section  3.4.1.  Two  subcases  were  considered: 
one  using  the  original  approximated  despun  platform  area,  and  the  second  using  varying  areas  (of 
five  and  ten  times  the  original  area)  to  illustrate  the  effects  of  having  a  larger  despun  platform  area 
on  the  semimajor  axis  of  the  satellite’s  orbit. 
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4-4- 1.1  Original  Despun  Platform  Area.  To  facilitate  comparison  between  the  purely 
specular  reflection  results  and  the  baseline  results,  the  respective  results  for  the  behavior  of  the 
semimajor  axis  over  time  are  overlaid  in  Figures  4.23-4.24. 

For  one  orbit,  as  shown  in  Figure  4.23,  the  simple  controller  yields  a  raise  in  semimajor  axis 
of  0.32216  km.  The  specular  reflection  case,  therefore,  exhibits  the  most  change  in  senrinrajor  axis 
thus  far,  and  one  would  expect  the  remaining  results  in  Section  4.4.1  to  corroborate  this.  Over  100 
orbits  A  a  is  4.0176  km,  and  over  1000  orbits  Aa  is  40.0064  km.  These  values  do  in  fact  exceed  the 
raises  in  semimajor  axis  presented  in  the  previous  sections  (for  the  simple  controller  and  nominal 
area) . 


x  1q4  a  versus  time  over  1  orbital  periods;  A  a  =0.32216  km 


Figure  4.23:  Simple  Controller  Case  Behavior  of  a  for  all  Specular  Reflection,  1  Orbit  and  an 

Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


44. 1.2  Varied  Despun  Platform  Area.  In  Section  4. 4. 1.1,  the  results  presented  were 
the  highest  raise  in  semimajor  axis  thus  far,  so  one  would  expect  the  trend  to  continue  for  larger 
despun  platform  areas.  This  section  addresses  the  results  when  area  factors  of  five  and  ten  were 
applied  to  the  despun  platform  area. 
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Figure  4.24:  Simple  Controller  Case  Behavior  of  a  for  all  Specular  Reflection,  1000  Orbits  and 

an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

For  an  area  factor  of  five,  the  results  in  Figure  4.25  exhibit  a  marked  improvement  over 
the  original  area  results  in  Section  4. 4. 1.1,  and  have  nearly  double  the  semimajor  axis  rise  of  the 
corresponding  all  absorption  case.  For  one  orbit,  Aa  is  0.5537  km.  For  100  orbits,  the  raise 
is  20.0511  bn  (the  all  absorption  case  was  10.4542  km,  by  comparison);  for  1000  orbits,  it  is 
200.5718  km  (the  all  absorption  case  was  104.3709  km). 

Figure  4.26  (1000  orbits)  was  generated  with  the  area  factor  set  equal  to  ten.  As  expected, 
the  raise  in  semimajor  axis  is  greater  in  this  case  than  for  both  the  nominal  area  and  with  an  area 
factor  of  five  for  the  same  scenario.  The  raise  in  semimajor  axis  over  one  orbit  is  0.84314  km,  which 
is  approximately  0.3  km  more  than  for  an  area  factor  of  five  and  0.5  km  more  than  for  the  nominal 
area.  This  result  is  also  the  best  for  any  one-orbit  case  so  far,  although  one  would  expect  the 
complex  case  to  generate  a  larger  change  in  a  since  the  controller  was  designed  for  purely  specular 
reflection. 

Over  100  orbits,  the  simple  controller  generated  a  Aa  of  40.0860  km,  double  that  of  the  area 
factor  of  five  and  ten  times  that  of  the  nominal  area  for  this  case.  For  1000  orbits,  the  Aa  is 
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Figure  4.25:  Simple  Controller  Case  Behavior  of  a  for  all  Specular  Reflection,  1000  Orbits  and 

an  Area  Factor  of  5,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

402.4059  km,  which  again  is  roughly  twice  that  of  the  area  factor  of  five  and  ten  times  that  of  the 
nominal  area  for  this  case.  Also,  as  for  the  one-orbit  case  above,  the  values  of  Aa  for  each  set  of 
orbits  are  the  highest  yet. 

4-4-%  Complex  Controller  Case.  This  section  covers  the  test  case  results  for  the  complex 
controller.  Again,  the  scenarios  used  all  specular  reflection.  Two  subcases  were  considered:  one 
with  the  nominal  despun  platform  area  and  the  second  with  area  factors  of  five  and  ten  applied  to 
the  nominal  area.  The  percent  improvement  over  the  simple  controller  is  approximately  42%. 

4-4-2.1  Original  Despun  Platform  Area.  The  satellite’s  original  effective  despun 
platform  area  was  used  to  compute  the  change  in  semimajor  axis  over  time,  as  illustrated  in 
Figure  4.27  for  1000  orbits.  Aa  for  one  orbit  is  0.34002  km,  which  is  a  slight  improvement  over  the 
simple  controller  case  in  Section  4. 4. 1.1.  Over  100  orbits,  the  raise  in  semimajor  axis  is  5.6209  km, 
and  for  1000  orbits  it  is  56.8227  km.  As  expected,  these  yield  the  best  results  so  far  for  the 
respective  number  of  orbits. 
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Figure  4.26:  Simple  Controller  Case  Behavior  of  a  for  all  Specular  Reflection,  1000  Orbits  and 

an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


Figure  4.27:  Complex  Controller  Case  Behavior  of  a  for  all  Specular  Reflection,  1000  Orbits  and 
an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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f.f.2,.2  Varied  Despun  Platform  Area.  This  section  addresses  the  changes  in  semi- 
major  axis  over  time  using  area  factors  of  five  and  ten  applied  to  the  nominal  despun  platform  area. 
For  an  area  factor  of  five,  the  semimajor  axis  exhibits  a  raise  of  0.6430  km  for  one  orbit,  which  is 
0.09  km  greater  than  for  the  simple  controller.  A  a  for  100  orbits  is  28.0751  krn  (about  8  km  more 
than  for  the  simple  controller)  and  for  1000  orbits  it  is  285.0813  km .  (roughly  85  km  more  than  for 
the  simple  controller). 

Figure  4.28  reveals  the  results  generated  when  an  area  factor  of  ten  is  applied  to  the  original 
despun  platform  area.  For  one  orbit,  the  change  in  semimajor  axis  is  1.0217  km,  which  is  about 
0.18  km  more  than  the  simple  controller  case  and  0.68  km  more  than  the  baseline  case.  Aa  for 
100  orbits  is  56.1473  km,  which  is  roughly  twice  that  of  the  complex  controller  for  an  area  factor 
of  five  and  1.5  times  the  simple  controller  of  the  same  case.  For  1000  orbits,  Aa  is  573.5517  km. 
As  expected,  the  results  for  1000  orbits  in  this  section  represent  the  greatest  change  in  semimajor 
axis  thus  far  (and  will  prove,  looking  ahead,  to  be  the  largest  Aa  of  all  the  cases). 


Figure  4.28:  Complex  Controller  Case  Behavior  of  a  for  all  Specular  Reflection,  1000  Orbits  and 
an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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4-5  Diffuse  Reflection  Case 


The  cases  in  this  section  implement  the  scenarios  for  all  diffuse  reflection;  the  diffuse  reflection 
coefficient  is  set  equal  to  one  and  the  remaining  coefficients  are  set  equal  to  zero.  As  in  the  previous 
sections,  the  performance  of  the  satellite  will  be  assessed  under  two  controllers  and  three  distinct 
area  options. 

4-5.1  Simple  Controller  Case.  The  results  for  the  simple  controller,  described  in  Sec¬ 
tion  3.4.1,  appear  in  this  section.  For  purposes  of  comparison,  the  results  are  again  overlaid  on  the 
baseline  case. 

4 .5. 1.1  Original  Despun  Platform  Area.  The  first  test  case  for  the  simple  controller 
employs  the  original  despun  platform  area.  For  one  orbit  (Figure  4.29),  the  raise  in  senrinrajor 
axis  is  0.31328  km .,  which  is  more  than  for  the  nominal  and  all  absorption  cases  but  less  than  for 
the  all  specular  reflection  case  (the  remaining  results  for  Aa  in  Section  4.5  will  follow  this  pattern 
of  being  slightly  less  than  the  specular  case).  A  a  for  100  orbits  is  3.4368  km  and  for  1000  orbits 
(Figure  4.30)  it  is  34.1497  km.. 

4-5. 1.2  Varied  Despun  Platform  Area.  To  illustrate  the  change  in  behavior  of  the 
semimajor  axis  over  time  for  different  despun  platform  areas,  this  section  applies  area  factors  of 
five  and  ten,  respectively,  to  the  nominal  area.  The  first  plot,  Figure  4.31,  depicts  the  semimajor 
axis  changes  over  1000  orbits  for  an  area  factor  of  five. 

For  one  orbit,  the  raise  in  semimajor  axis  is  0.50933  km.,  for  100  orbits  it  is  17.1353  km  and 
for  1000  orbits  it  is  171.3268  km..  As  expected,  these  values  lie  (respectively  per  number  of  orbits) 
between  the  nominal  case  and  the  specular  reflection  case  for  an  area  factor  of  five. 

An  area  factor  of  ten  yields  better  results;  the  raise  in  semimajor  axis  for  one  orbit  is 
0.75438  km,  for  100  orbits  it  is  34.2684  km  and  for  1000  orbits  it  is  343.6828  km  (as  illustrated 
in  Figure  4.32).  For  1000  orbits,  the  diffuse  reflection  case  yields  a  slight  improvement  over  the 
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Figure  4.29:  Simple  Controller  Case  Behavior  of  a  for  all  Diffuse  Reflection,  1  Orbit  and  an  Area 

Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


1  q«  a  versus  time  over  1 000  orbital  periods;  A  a  =34.1 497  km 


Figure  4.30:  Simple  Controller  Case  Behavior  of  a  for  all  Diffuse  Reflection,  1000  Orbits  and  an 

Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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Figure  4.31:  Simple  Controller  Case  Behavior  of  a  for  all  Diffuse  Reflection,  1000  Orbits  and  an 

Area  Factor  of  5,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

nominal  case  (about  25  km),  a  marked  improvement  over  the  all  absorption  case  (roughly  134  km) 
and  a  small  disadvantage  when  compared  to  the  specular  reflection  case  (by  about  59  km). 

4-5.2  Complex  Controller  Case.  This  section  addresses  the  changes  in  semimajor  axis 
for  the  diffuse  case  when  the  complex  controller  is  used  versus  the  simple  controller.  Again,  two 
subcases  were  considered:  one  with  the  original  estimated  despun  platform  area  and  the  second 
with  area  factors  of  five  and  ten  applied  to  the  original  area.  The  percent  difference  from  using  the 
complex  controller  over  the  simple  controller  is  —1.7%. 

4- 5. 2.1  Original  Despun  Platform  Area.  The  nominal  area  results  are  illustrated 
for  1000  orbits  in  Figure  4.33.  Similarly  to  the  all  absorption  case,  the  resulting  increases  for  the 
semimajor  axis  using  the  complex  controller  are  slightly  less  than  using  the  simple  controller. 

The  raise  in  semimajor  axis  for  one  orbit  is  0.31559  km ,  which  is  0.0023  km.  more  than  for 
the  simple  controller.  For  100  orbits  Aa  is  3.3594  km  and  for  1000  orbits  Aa  is  33.5575  km.  For  all 
diffuse  reflection,  the  complex  controller  cases  yield  a  slightly  lower  A  a  than  their  corresponding 
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Figure  4.32:  Simple  Controller  Case  Behavior  of  a  for  all  Diffuse  Reflection,  1000  Orbits  and  an 

Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 

simple  controller  cases  for  100  and  1000  orbits  (note,  the  differences  between  the  controllers’  results 
for  the  diffuse  case  are  minimal  as  compared  to  the  all  absorption  case). 

4- 5.2.2  Varied  Despun  Platform  Area.  The  results  presented  in  this  section  use  two 
different  area  multiplication  factors  (five  and  ten)  applied  to  the  original  despun  platform  area. 

For  an  area  factor  of  five,  the  resulting  changes  in  semimajor  axis  are  less  than  for  the 
corresponding  simple  controller  case.  For  one  orbit,  Aa  is  0.52087  /cm;  for  100  orbits,  Aa  is 
16.7510  km;  for  1000  orbits,  Aa  is  168.2071  km. 

Figure  4.34  illustratea  the  change  in  semimajor  axis  versus  time  for  an  area  factor  of  ten 
applied  to  the  nominal  despun  platform  area.  While  the  results  still  vary  dramatically  from  the 
baseline  results,  Aa  is  slightly  less  in  this  case  than  for  the  corresponding  simple  controller  case. 
For  one  orbit,  Aa  is  0.77747  km.  For  100  orbits,  Aa  is  33.4959  km  and  for  1000  orbits,  Aa  is 
337.7322  km. 
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Figure  4.33:  Complex  Controller  Case  Behavior  of  a  for  all  Diffuse  Reflection,  1000  Orbits  and 

an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


Figure  4.34:  Complex  Controller  Case  Behavior  of  a  for  all  Diffuse  Reflection,  1000  Orbits  and 

an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 
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4-6  Specular  and  Diffuse  Reflection  Case 


This  section  addresses  the  case  in  which  the  satellite’s  diffuse  and  specular  coefficients  of 
reflection  are  each  set  to  equal  one-half  and  the  coefficient  of  absorption  is  set  to  zero.  One  would 
expect  the  results  for  this  case  (hereafter  referred  to  as  the  mixed  reflection  case)  to  fall  between 
the  results  from  the  previous  two  sections  (purely  diffuse  and  purely  specular  reflections).  The 
resultant  force  vector  for  the  mixed  case  will  lie  between  the  resultant  force  vectors  for  the  purely 
specular  and  purely  diffuse  cases,  which  leads  to  the  change  in  senrinrajor  axis  being  in  between 
the  two  cases  as  well. 

4-6.1  Simple  Controller  Case.  The  simple  controller,  described  in  Section  3.4.1,  is  imple¬ 
mented  in  this  section.  The  mixed  reflection  case  of  half  specular  and  half  diffuse  reflection  was 
used  to  produce  the  plots  in  the  following  sections.  Two  subcases  were  considered:  one  using  the 
original  despun  platform  area,  and  the  second  using  varying  areas  to  illustrate  the  effects  of  larger 
despun  platform  areas  on  the  semimajor  axis. 

4-6. 1.1  Original  Despun  Platform  Area.  The  results  in  this  section  use  the  original 
despun  platform  area.  To  facilitate  straightforward  comparison,  the  results  from  this  section  were 
plotted  on  the  same  graphs  as  the  baseline  results  from  Section  4.2.1.  Figures  4.35-4.36  illustrate 
the  original  area  results  for  the  mixed  reflection  case. 

For  one  orbit,  the  raise  in  semimajor  axis  is  0.31772  km. ,  which  does  lie  between  the  corre¬ 
sponding  purely  specular  (0.32216  km)  and  purely  diffuse  (0.31328  km)  values  as  expected.  For 
100  orbits,  A  a  is  3.7267  km ;  for  1000  orbits,  it  is  37.0899  km.  These  values  also  lie  between  the 
corresponding  purely  specular  and  diffuse  reflection  results. 

4-6. 1.2  Varied  Despun  Platform  Area.  The  results  in  this  section  use  the  two 
area  factors  applied  to  the  original  despun  platform  area.  The  first  area  factor,  five,  is  used  in 
Figure  4.37. 
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1  q*  a  versus  time  over  1  orbital  periods;  A  a  =0.31  111  km 


Figure  4.35:  Simple  Controller  Case  Behavior  of  a  for  50%  Diffuse  and  50%  Specular  Reflection, 

1  Orbit  and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of  a. 


Figure  4.36:  Simple  Controller  Case  Behavior  of  a  for  50%  Diffuse  and  50%  Specular  Reflection, 

1000  Orbits  and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of 
a. 
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For  one  orbit,  the  resulting  change  in  semimajor  axis  is  0.53151  km.  Over  100  orbits,  the 


total  raise  is  18.5895  km  and  for  1000  orbits  it  is  185.9499  km.  Again,  these  values  fall  between 
the  purely  specular  and  diffuse  reflection  values  for  the  same  case. 


Figure  4.37:  Simple  Controller  Case  Behavior  of  a  for  50%  Diffuse  and  50%  Specular  Reflection, 

1000  Orbits  and  an  Area  Factor  of  5,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of 
a. 


Figure  4.38  illustrates  the  resulting  changes  in  semimajor  axis  when  an  area  factor  of  ten  is 
applied  to  the  original  despun  platform  area.  For  one  orbit,  Aa  is  0.79876  km.  Over  100  orbits, 
it  is  37.1804  km  and  over  1000  orbits  it  is  373.0811  km.  Once  again,  as  expected,  these  values  for 
the  raise  in  semimajor  axis  fall  between  the  purely  specular  and  purely  diffuse  results. 

4-6.2  Complex  Controller  Case.  This  section  presents  the  results  for  the  complex  con¬ 
troller  described  in  Section  3.4.2.  Again,  two  subcases  were  considered:  one  using  the  original 
despun  platform  area  and  the  second  using  two  area  factors  applied  to  the  original  despun  plat¬ 
form  area.  The  percent  difference  from  using  the  complex  controller  over  the  simple  controller  is 
approximately  22%. 


84 


Figure  4.38:  Simple  Controller  Case  Behavior  of  a  for  50%  Diffuse  and  50%  Specular  Reflection, 

1000  Orbits  and  an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of 
a. 

4-6.2. 1  Original  Despun  Platform  Area.  Figure  4.39  illustrates  the  resulting  changes 
in  semimajor  axis  using  the  complex  controller  and  the  original  despun  platform  area  of  the  satellite. 

As  expected,  the  values  of  Aa  for  this  particular  case  lie  between  the  purely  specular  and 
purely  diffuse  reflection  for  the  corresponding  number  of  orbits.  For  this  scenario,  over  one  orbit, 
the  raise  in  semimajor  axis  is  0.3278  km.  Over  100  orbits,  Aa  is  4.4885  km  and  for  1000  orbits  it 
is  45.1932  km.  (approximately  10  km  more  than  the  simple  controller  case  for  mixed  reflection  and 
original  area). 

4-6. 2. 2  Varied  Despun  Platform  Area.  This  section  details  the  results  when  area 
factors  of  five  and  ten  are  applied  to  the  original  despun  platform  area. 

For  an  area  factor  of  five  over  one  orbit,  the  raise  in  semimajor  axis  is  0.58193  km.  A  a  is 
22.4093  km.  for  100  orbits  and  226.6733  km  for  1000  orbits.  These  values  again  fall  between  the 
corresponding  raises  in  semimajor  axis  for  the  specular  and  diffuse  cases. 
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Figure  4.39:  Complex  Controller  Case  Behavior  of  a  for  50%  Diffuse  and  50%  Specular  Reflection, 
1000  Orbits  and  an  Area  Factor  of  1,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of 
a. 

Figure  4.40  shows  the  results  for  an  area  factor  of  ten  applied  to  the  original  despun  platform 
area.  For  one  orbit,  the  change  in  semimajor  axis  is  0.8996  km.  A  a  is  44.8291  km .  for  100  orbits 
and  is  455.0825  km.  for  1000  orbits.  These  also  fall  between  the  values  for  the  specular  and  diffuse 
reflection  cases. 


4-7  Summary  of  Results 

A  tabular  summary  of  the  results  from  the  thesis  is  presented  in  Table  4.1.  In  Table  4.1,  note 
that  for  all  cases,  purely  specular  reflection  results  in  the  largest  change  in  semimajor  axis.  For  the 
simple  controller,  the  purely  diffuse  case  produces  the  third  best  results  for  semimajor  axis  change, 
primarily  because  the  half-diffuse  and  half-specular  mixed  case  results  lie  in  between  the  purely 
specular  and  purely  diffuse  results.  Notably,  the  complex  controller  yields  less  desirable  results 
than  the  simple  controller  for  the  purely  diffuse  case,  which  is  likely  due  to  the  complex  controller 
being  designed  to  take  advantage  of  completely  specular  reflection.  In  fact,  the  nominal  case  for 
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Figure  4.40:  Complex  Controller  Case  Behavior  of  a  for  50%  Diffuse  and  50%  Specular  Reflection, 
1000  Orbits  and  an  Area  Factor  of  10,  Overlaid  with  the  Corresponding  Baseline  Case  Behavior  of 
a. 

the  complex  controller  yielded  a  longer  rise  in  semimajor  axis  than  the  diffuse  case.  The  satellite’s 
nominal  reflection  case,  with  each  coefficient  set  to  one-third,  yields  results  that  are  slightly  less 
than  the  purely  diffuse  case  for  the  simple  controller.  Finally,  as  expected,  the  all  absorption  case 
generates  the  least  raise  in  semimajor  axis  over  all  the  cases. 

Overall,  then,  the  order  of  effectiveness  is: 

1.  Complex  controller  with  crs  =  1 

2.  Complex  controller  with  =  0.5  and  crs  =  0.5 

3.  Simple  controller  with  crs  =  1 

4.  Simple  controller  with  cr<i  =  0.5  and  crs  =  0.5 

5.  Complex  controller  with  CRd  =  crs  =  CRa  =  5 

6.  Simple  controller  with  CRd  =  1 

7.  Complex  controller  with  CRd  =  1 

8.  Simple  controller  with  CRd  =  crs  =  crq  =  | 

9.  Simple  controller  with  CRa  =  1 

10.  Complex  controller  with  end  =  1 
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Table  4.1:  Summary  of  Results. 


Case 

Aa  (km),  the  change  in  semimajor  axis  from  the  epoch  value 

1  Orbit 

100  Orbits 

1000  Orbits 

Area  Factor 

1 

5 

10 

1 

5 

10 

1 

5 

10 

Baseline 

0.3546 

0.3546 

0.3546 

0.0129 

0.0129 

0.0129 

0.0105 

0.0105 

0.0105 

Nominal 

Simple 

Complex 

0.3098 

0.4920 

0.7197 

3.1860 

15.8804 

31.7532 

31.6495 

158.7181 

318.2681 

0.3147 

0.5163 

0.7684 

3.4266 

17.0842 

34.1622 

34.3087 

172.1240 

345.1985 

Absorption 

Simple 

Complex 

0.2940 

0.4129 

0.5614 

2.1010 

10.4542 

20.9048 

20.7914 

104.3709 

209.1176 

0.2884 

0.3851 

0.5056 

1.2974 

6.4357 

12.8627 

12.6168 

63.2965 

126.7702 

Specular 

Simple 

Complex 

0.3222 

0.5537 

0.8431 

4.0176 

20.0511 

40.0860 

40.0064 

200.5718 

402.4059 

0.3400 

0.6430 

1.0217 

5.6209 

28.0751 

56.1473 

56.8227 

285.0813 

573.5517 

Diffuse 

Simple 

Complex 

0.3133 

0.5093 

0.7544 

3.4368 

17.1353 

34.2684 

34.1497 

171.3268 

343.6828 

0.3156 

0.5209 

0.7775 

3.3594 

16.7510 

33.4959 

33.5575 

168.2071 

337.7322 

Mix 

Simple 

Complex 

0.3177 

0.5315 

0.7988 

3.7267 

18.5895 

37.1804 

37.0899 

185.9499 

373.0811 

0.3278 

0.5819 

0.8996 

4.4885 

22.4093 

44.8291 

45.1932 

226.6733 

455.0825 

As  the  main  objective  of  the  thesis  is  to  raise  the  semimajor  axis  of  the  DSCS  II  satellite 
into  a  disposal  orbit,  it  is  beneficial  to  address  the  length  of  time  required  to  do  this.  Based  on  a 
linear  extrapolation  of  the  nominal  absorption  and  reflection  results  for  the  complex  controller  and 
an  area  factor  of  one  (so  as  to  use  the  satellite’s  actual  estimated  area),  the  time  required  to  raise 
the  DSCS  II  400  km.  would  be  about  33  years,  or  about  12000  orbits. 

Comparatively,  the  time  required  is  roughly  divisible  by  the  area  multiplication  factor  used 
(keeping  a  constant  mass  of  611  kg);  a  change  of  400  km  could  be  reached  in  approximately  6.6  years 
if  the  area  of  the  satellite  was  five  times  the  original  area,  and  only  3.3  years  if  it  was  ten  times 
more.  Hence,  perhaps  the  best  takeaway  from  the  results  chapter  is  that  even  a  small  change  in 
apparent  area  (area  of  the  satellite  facing  the  Sun)  yields  a  significantly  shorter  time  required  to 
raise  the  satellite  into  a  disposal  orbit.  This  point  is  reinforced  in  Figure  4.41,  which  illustrates  the 
effects  of  changing  areas  and  surface  reflection  properties  on  the  time-to-disposal.  Incidently,  this 
effect  is  not  a  surprise,  as  the  control  force  available  is  linearly  proportional  to  the  controlled  area. 


Time  to  400  km  orbit  raise  for  varying  methods 


Figure  4.41:  Time  (in  Years)  to  400  km  versus  Case  for  Complex  Controller.  The  plot  illustrates 
how  much  time  is  required  to  raise  the  orbit  for  differing  area  factors  and  coefficients  of  absorption 
and  reflection.  The  case  type  and  area  factor  are  labeled  on  the  x-axis. 
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The  complex  controller  case  with  an  area  factor  of  ten  and  purely  specular  reflection  yields 
the  fastest  time-to-disposal  (approximately  two  years).  For  the  original  area,  the  specular  case 
again  provides  the  least  amount  of  time  required  for  disposal,  which  reflects  the  optimization  of  the 
complex  controller  for  pure  specular  reflection.  Hence,  the  nominal  (“best  guess”  for  the  satellite) 
cases  yield  more  desirable  results  than  both  the  pure  absorption  and  diffuse  reflection  cases.  These 
results  are  also  consistent  with  the  equations  for  the  perturbed  acceleration  due  to  SRP  from 
Chapter  III,  since  pure  specular  reflection  results  in  the  highest  potential  force  per  unit  area  (as 
illustrated  in  Figures  3.3  and  3.4). 

If  instead  the  satellite  is  best  modeled  as  a  mix  of  specular  and  diffuse  reflection,  which  is 
most  likely  a  relatively  good  approximation  also  as  the  majority  of  the  satellite  is  composed  of 
highly  polished  surfaces  (the  solar  arrays),  the  time-to-disposal  decreases  from  the  nominal  case 
by  about  8  years.  If  the  area  of  the  satellite’s  despun  platform  was  increased  by  a  factor  of  five, 
from  roughly  5  in 2  to  25  m2  (or  2.7  m  by  6.1  m  versus  1.2  in  by  2.7  in),  both  the  nominal  and 
mixed  reflection  cases  would  enable  the  satellite  to  reach  a  400  km  raise  in  approximately  5  or  6 
years.  Although  a  retrofit  to  the  in-orbit  DSCS  II  F-13  modeled  in  this  thesis  would  be  costly  (and 
unnecessary,  as  if  astronauts  were  to  retrofit  it,  they  might  as  well  dismantle  it  and  bring  it  down 
from  orbit),  this  is  an  interesting  consideration  for  future  satellites.  Instead  of  thrusters  and  fuel, 
future  satellites  could  store  a  small  deployable  solar  sail  with  which  to  maneuver  (using  electric 
power)  into  a  disposal  orbit  at  the  end  of  their  mission  lives. 

In  short,  the  two  factors  which  have  the  most  effect  on  the  amount  of  time  it  takes  to  raise  a 
satellite  into  a  disposal  orbit  are  the  area-to-mass  ratio  and  the  optical  properties  of  the  satellite 
surfaces.  There  does  not  appear  to  be  a  “preferred”  area,  as  the  disposal  time  versus  the  area  was 
essentially  linear  throughout  the  region  simulated.  Additional  comments  on  the  results,  as  well  as 
suggested  further  work,  can  be  found  in  Chapter  V. 
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V.  Conclusions  and  Recommendations 


5.1  Summary  and  Conclusions 

This  thesis  has  demonstrated  that  solar  radiation  pressure,  as  the  most  significant  non- 
gravitational  perturbation,  does  have  the  potential  to  raise  geosynchronous  satellite  into  disposal 
orbits.  Using  a  controller,  the  satellite  can  be  maneuvered  to  optimize  the  direction  of  incident 
solar  radiation  on  its  surface,  thus  using  the  force  from  SRP  to  raise  the  satellite  into  a  disposal  or¬ 
bit.  Several  factors  control  the  accuracy  of  the  model;  conical  versus  cylindrical  Earth  shadow  (the 
satellite  spends  on  the  order  of  5%  and  1%  of  the  orbital  period  in  shadow,  respectively  [19,20]), 
constant  versus  varying  solar  flux  (about  a  7%  difference  between  minimum  and  maximum  val¬ 
ues),  constant  versus  changing  cross-sectional  area  of  the  satellite,  and  the  optical  properties  of  the 
satellite  surface  (absorption,  reflection,  and  transmission).  A  derivation  of  the  SRP  model  incor¬ 
porating  a  cylindrical  Earth  shadow  and  time-varying  solar  flux  was  developed  in  Chapter  III.  In 
turn,  numerical  and  analysis  results  were  presented  in  Chapter  IV. 

Despite  the  potential  of  SRP,  it  is  clear  that  it  may  not  be  a  practical  solution  for  existing 
spacecraft.  Under  the  best-performing  surface  properties  (purely  specular  reflection)  and  using  the 
best  controller,  the  time-to-disposal  is  on  the  order  of  20  years  for  the  DSCS  II  satellite  considered 
in  this  work.  The  satellite,  however,  most  likely  retains  surface  characteristics  similar  to  either  the 
nominal  (all  coefficients  of  reflection  and  absorption  set  equal  to  1/3)  or  mixed  (the  coefficients  of 
specular  and  diffuse  reflection  each  set  equal  to  1/2)  cases.  As  the  complex  controller  provides  the 
largest  change  in  semimajor  axis  for  these  two  cases,  reference  Figure  4.41  for  the  corresponding 
times-to-disposal  of  approximately  33  years  (nominal)  and  25  years  (mixed).  For  larger  area-to- 
mass  ratios,  the  time  required  to  raise  the  DSCS  II  F-13  satellite  400  km  is  lowered  as  one  over  the 
area  factor.  Hence,  a  small  change  in  area  has  a  large  effect  on  the  time-to-disposal  of  the  satellite, 
as  does  small  changes  in  material  properties  of  the  satellite’s  surfaces.  This  offers  hope  that  new 
systems  designed  to  take  advantage  of  SRP  could  achieve  much  better  disposal  times. 
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5.2  Recommendations  for  Future  Work 


Future  research  related  to  this  thesis  would  prove  to  be  beneficial  in  estimating  to  a  higher 
order  of  accuracy  the  time  required  to  move  the  DSCS  II  satellites  into  disposal  orbits.  Several 
simplifications  were  made  in  the  thesis  which,  if  included  in  later  work,  would  yield  a  more  likely 
estimate. 

A  conical  shadow  model,  while  more  difficult  to  implement,  is  a  more  accurate  representation 
of  the  Earth’s  shadow  than  the  cylindrical  model.  Modeling  the  shadow  as  such  would  enable  the 
researcher  to  study  the  transition  effects  of  the  spacecraft  through  the  regions  of  the  umbra  and 
penumbra.  The  use  of  a  cylindrical  model  neglects  these  transition  regions;  it  models  the  shadow  as 
having  only  the  dark,  or  umbra,  regions.  This,  in  turn,  causes  the  solar  radiation  pressure  to  only 
appear  “off”  (satellite  in  the  umbra)  or  “on”  (satellite  in  sunlight).  The  amount  of  time  spacecraft 
spend  out  of  the  shadow  region  directly  correlates  to  the  amount  of  time  needed  to  raise  its  orbit, 
and  a  better  model  for  the  shadow  would  produce  a  more  accurate  value  for  this  time-to-disposal. 
An  alternate  approach  would  be  to  model  the  cylindrical  shadow  with  a  reduced  size;  this  would 
test  the  importance  of  the  shadow  as  it  affects  the  solution. 

Also  of  importance  is  the  method  developed  by  Coverstone  et.  al.  for  the  complex  controller. 
In  additions  to  the  modifications  made  in  this  thesis  (reference  Chapter  I)  to  account  for  the 
inclination  of  the  Earth  with  respect  to  the  ecliptic  plane,  future  work  could  further  modify  her 
model  to  include  absorption,  diffuse  reflection  and  transmission  of  the  incident  solar  ray  on  the 
surface  of  the  satellite.  As  her  model  only  includes  specular  reflection,  the  results  are  certainly 
skewed  to  favor  the  case  with  purely  specular  reflection.  The  inclusion  of  more  complex  surface 
properties  of  the  satellite  would  generate  potentially  different  results. 

If  continual  interest  in  raising  the  orbits  of  the  DSCS  II  satellites  were  present,  future  work 
should  include  better  data  (if  possible)  on  the  properties  of  the  satellites.  Although  the  properties 
of  the  satellite,  to  include  mass,  area,  dimensions  of  the  parabolic  antennae  and  coefficients  of 
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reflection  and  absorption,  were  reasonably  know  at  the  time  of  launch  more  than  20  years  ago,  the 
current  values  are  not  as  readily  available.  As  the  results  show,  these  values  can  have  a  profound 
impact  on  the  time  required  to  dispose  of  the  satellite  into  a  higher  orbit.  Of  particular  interest 
is  the  parabolic  antennae;  this  thesis  uses  a  flat  plate  approximation  for  them,  and  having  the 
appropriate  dimensions  available  would  enable  a  more  realistic  propagation. 

Applying  this  method  to  other  aging  satellites  in  orbit  might  also  be  interesting,  particularly 
for  satellites  of  different  area-to-mass  ratios  and  those  in  different  orbits  than  geosynchronous.  An 
important  system-level  trade-off  consideration  would  be  the  fuel  mass  needed  for  disposal  versus  the 
mass  required  to  store,  deploy  and  control  a  large  area  “disposal  sail” .  Of  course,  if  used  for  satellites 
in  lower  orbits  (less  than  800  km ),  drag  will  quickly  become  the  highest-order  non-gravitational 
perturbation  instead  of  SRP. 

The  author  believes  that  additional  study  of  using  SRP  to  dispose  of  satellites  otherwise 
incapable  of  being  moved  into  disposal  orbits  will  remain  of  interest.  As  long  as  the  demand 
for  newer  (or  simply  more)  satellites  in  high-value  orbits  such  as  the  geosynchronous  belt  exists, 
researchers  will  continue  to  be  tasked  to  address  how  to  dispose  of  satellites  that  have  reached  the 
end  of  their  mission  life.  For  existing  dual-spin  satellites  like  the  DSCS  II  modeled  in  this  thesis, 
however,  this  method  is  unlikely  a  practical  option.  Instead,  future  work  could  be  performed  to 
investigate  applications  with  a  three-axis  satellite  with  large  solar  array  panels  (much  like  the 
current  constellation  of  DSCS  III  satellites).  The  larger  area  may  be  able  to  be  controlled  to  both 
provide  sufficient  power  and  perform  orbit  raising.  Another  potential  use  is  for  geosynchronous 
satellite  stationkeeping;  if  the  design  of  the  satellite  included  a  sail-type  device,  it  could  be  used 
for  both  stationkeeping  throughout  the  satellite’s  mission  life  and  for  disposal  at  end-of-life. 
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Appendix  A.  DSCS  II  TLE  and  Satellite  Parameters 


A.l  Interpreting  the  Two-Line  Element  Set 

The  most  common  method  for  communicating  orbit  information  is  the  TLE.  This  format, 
used  by  NASA  and  the  North  American  Aerospace  Defense  Command  (NORAD),  consists  of  three 
lines  of  information  that  can  be  deciphered  to  obtain  the  properties  of  the  satellite  and  its  current 
position  in  orbit  about  the  Earth. 

A  TLE  has  one  line  consisting  of  a  24-character  name  (line  0)  and  two  69-character  lines 
of  data  (lines  1  and  2).  The  only  valid  characters  in  the  lines  are  the  numbers  0-9,  the  capital 
letters  A-Z,  the  period,  the  space,  and  the  plus  and  minus  signs.  There  are  further  restrictions  on 
the  elements  of  each  line,  which  are  detailed  on  the  CelesTrak  website  ( http://celestrak.com ).  The 
generic  format  of  a  TLE  is 

AAAAAAAAAAAAAAAAAAAAAAAA 

1  NNNNNU  NNNNNAAA  NNNNN . NNNNNNNN  + . NNNNNNNN  +NNNNN-N  +NNNNN-N  N  NNNNN 

2  NNNNN  NNN . NNNN  NNN . NNNN  NNNNNNN  NNN . NNNN  NNN . NNNN  NN . NNNNNNNNNNNNNN 

Tables  A.l  and  A. 2  define  each  of  the  individual  Helds  for  lines  1  and  2,  respectively  [15,21]. 

The  code  presented  in  Appendix  C  requires  the  input  of  most  of  the  parameters  in  the  TLE, 
and  uses  them  to  calculate  the  classical  orbital  elements  at  epoch.  The  TLE  used  for  this  thesis 
is  [28] 


OPS  9443 

1  11621U  T9098A  05202.40912884  -.00000117  00000-0  10000-3  0  8408 

2  11621  014.4997  018.8800  0005671  321.5388  038.4719  00.95598989542422 
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Table  A.l:  TLE  Line  1  Description. 


Column 

Description 

01 

Line  number  of  element  data 

03-07 

Satellite  number  (from  NORAD  number) 

08 

Classification  of  TLE  (U=Unclassified  and  C=Classified) 

10-11 

International  designator  (last  two  digits  of  launch  year) 

12-14 

International  designator  (launch  number  of  the  year) 

15-17 

International  designator  (piece  of  the  launch) 

19-20 

Epoch  year  (only  the  last  two  digits) 

21-32 

Epoch  (day  of  the  year  and  fraction  of  the  day) 

34-43 

First  time  derivative  of  the  mean  motion  or  ballistic  coefficient* 

45-52 

Second  time  derivative  of  the  mean  motion  (decimal  point  assumed) 

54-61 

B*  drag  term  (decimal  point  assumed)  or  radiation  pressure  coefficient^ 

63 

Ephemeris  type 

65-68 

Element  number 

69 

Checksum  value 

*depends  on  ephemeris  type 
1 B *  if  GP4  general  perturbation  theory  was  used 


Table  A. 2:  TLE  Line  2  Description. 


Column 

Description 

01 

Line  number  of  element  data 

03-07 

Satellite  number 

09-16 

Inclination  (degrees) 

18-25 

Right  ascension  of  the  ascending  node  (degrees) 

27-33 

Eccentricity  (decimal  point  assumed) 

35-42 

Argument  of  perigee  (degrees) 

44-51 

Mean  anomaly  (degrees) 

53-63 

Mean  motion  (revolutions  per  day) 

64-68 

Revolution  number  at  epoch  (revolutions) 

69 

Checksum  value 
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A. 2  DSCS  II F-13  Data 


The  DSCS  satellite  system,  of  which  a  typical  Phase  II  satellite  is  illustrated  in  Figure  A.l,  is 
a  United  States  military  satellite  constellation  that  was  placed  in  a  geosynchronous  orbit  to  provide 
secure  voice  and  data  communications  to  the  military  user.  The  DSCS  system  succeeded  the  Initial 
Defense  Satellite  Communications  System  (IDSCS)  program,  which  launched  its  first  satellite  into 
orbit  in  1966. 


Several  reputable  sources  are  available  for  determining  the  characteristics  of  the  DSCS  satel¬ 
lites,  however,  they  often  do  not  agree  on  specifics.  Jane’s  All  the  World’s  Aircraft,  for  example, 
assumes  that  each  DSCS  satellite  has  the  same  specifications.  NASA  gives  details  for  each  satellite 
launched,  but  leaves  out  some  necessary  details.  The  specifications,  as  merged  from  both  sources 
for  the  DSCS  II  satellite  used  in  the  thesis,  are  outlined  in  Table  A. 3  [23,37]. 


Table  A. 3:  DSCS  II  F-13  Specifications. 


Parameter 

Value 

Units 

mass 

611 

kg 

Spin-Stabilized  Platform 

height 

1.8288 

m 

diameter 

2.7432 

m 

surface  area 

5.0168 

2 

m 

Despun  Platform 

height 

1.2192 

m 

width 

2.7432 

m 

surface  area' 

3.3445 

2 

m 

t  estimated  approximate  surface  area 


96 


mm 


jrrrrr^ 


Figure  A.l:  DSCS  II  Satellite  [1].  Photo  courtesy  of  NASA  JPL. 
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Appendix  B.  Aphelion  Almanac 


The  Earth  reaches  aphelion  around  4  July  annually.  The  variance  in  date  per  year  is  well  doc¬ 
umented  in  the  Astronomical  Almanac,  published  by  the  United  States  Naval  Observatory’s  As¬ 
tronomical  Applications  Department.  Table  B.l  is  taken  from  the  larger  database,  which  can  be 
found  on  their  website.  Here,  d,  h  and  m  indicate  day,  hour  and  minute,  respectively,  in  Universal 
Time  [31]. 


Table  B.l:  Earth  Aphelion  Dates  for  Years 
2005-2020. 


Year 

Month 

d 

h 

Day  of  Year 

2005 

July 

5 

05 

186 

2006 

July 

3 

23 

184 

2007 

July 

7 

00 

188 

2008 

July 

4 

08 

186 

2009 

July 

4 

02 

185 

2010 

July 

6 

11 

187 

2011 

July 

4 

15 

185 

2012 

July 

5 

03 

187 

2013 

July 

5 

15 

186 

2014 

July 

4 

00 

185 

2015 

July 

6 

19 

187 

2016 

July 

4 

16 

186 

2017 

July 

3 

20 

184 

2018 

July 

6 

17 

187 

2019 

July 

4 

22 

185 

2020 

July 

4 

12 

186 
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Appendix  C.  Matlab®  Source  Code 


The  Matlab®  code  is  presented  in  alphabetical  order  by  m-file.  To  facilitate  better  understanding 
of  the  overall  code  hierarchy,  reference  Figure  C.l. 


Equcode_wrapper.m 


juldatet.m 


~sr 

▼ 


earthsunt.m 


srpaccel.m 


Figure  C.l:  Basic  Matlab®  m-file  Structure. 
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C.l  CALCEA.m 


7.  Solves  for  eccentric  anomaly,  E  from  a  given  mean  anomaly,  M 
7.  and  eccentricty,  e.  Performs  a  simple  Newton-Raphson  iteration 

7. 

7,  M  must  be  in  radians  and  E  is  returned  in  radians. 

7. 

7.  Default  tolerances  is  10~-8  rad. 

7. 

7.  E  =  CalcEA(M,e)  uses  default  tolerances 

7. 

7,  E  =  CalcEA(M,e,tol)  will  use  a  user  specified  tolerance,  tol 

7. 

7.  Richard  Rieber 
7.  1/23/2005 

7.  E-mail  problems/suggestions  rrieber@gmail.com 

function  E  =  CalcEA(M,e,tol) 

7»Checking  for  user  inputed  tolerance 
if  nargin  ==  2 

7oUsing  default  value 
tol  =  10~-8; 
elseif  nargin  >  3 

error (’Too  many  inputs.  See  help  CalcE’) 
elseif  nargin  <  2 

error (’Too  few  inputs.  See  help  CalcE’) 

end 

Etemp  =  M;  ratio  =  1;  while  abs (ratio)  >  tol 
f_E  =  Etemp  -  e*sin(Etemp)  -  M; 
f_Eprime  =  1  -  e*cos (Etemp) ; 
ratio  =  f_E/f_Eprime; 
if  abs (ratio)  >  tol 

Etemp  =  Etemp  -  ratio; 

else 

E  =  Etemp; 

end 

end 


C.2  earthsun.m 

function  r_ES_vec  =  earthsun(epoch_yr ,epoch_day,epoch_frac) 

7,  This  function  calculates  the  position  vector  between  the 

7.  earth  and  the  sun  for  the  calculated  JD  in  UTC .  Note  that  this  algorithm 
7.  is  only  valid  until  the  year  2050  due  to  truncation  (see  Vallado  pg  265) . 


7.  note  all  of  these  values  must  be  in  degrees,  not  radians! ! ! 
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JD  =  juldate(epoch_yr,epoch_day,epoch_frac) ; 

7.  JD  is  in  ref  to  UTC  (let  JD_UTC  =  JD_UT1) 

%  find  the  number  of  Julian  centuries 

T_UT1  =  (JD  -  2451545. 0)/36525;  7,  Vallado  pg  188 

lambda_M_sundeg  =  280.460  +  36000 . 770*T_UT1 ;  7,  deg 

T_TDB  =  T_UT1 ; 

M_sundeg  =  357.5277233  +  35999 . 05034*T_TDB ;  1  deg 

lambda_ecliptic  =  lambda_M_sundeg  +  1 . 914666471*sind(M_sundeg)  + 

0 . 019994643*sind(2*M_sundeg) ;  7«  deg 

r_sun  =  1.000140612  -  0 . 016708617*cosd(M_sundeg)  -... 

0 . 000139589*cos  (2*M_sundeg)  ;  '/,  AU 

epsilon  =  23.439291  -  0 . 0130042*T_TDB ;  7.  deg 

7.  the  vector  referred  to  the  epoch  of  date  is  (in  AU) 
r_ES_vecAU  =  [r_sun*cosd(lambda_ecliptic) ; . . . 

r_sun*cosd(epsilon) *sind(lambda_ecliptic) ; . . . 
r_sun*sind (epsilon) *sind(lambda_ecliptic)]  ; 

r_ES_vec  =  149598000 . *r_ES_vecAU;  7o  Sun-Earth  vector  in  km 

7.  also  calculate  the  right  ascension  and  declination  of  the  sun  (can  check 
7.  with  Astronomical  Almanac) 

alpha_sun  =  atand(cosd(epsilon) *tand(lambda_ecliptic) ) ;  delta_sun 
=  asind(sind(epsilon)*sind(lambda_ecliptic)) ; 

return 


C.3  earthsunt.m 

function  r_ES_vec  =  earthsunt(t) 

7,  This  function  calculates  the  position  vector  between  the 

7.  earth  and  the  sun  for  the  calculated  JD  in  UTC.  Note  that  this  algorithm 
7.  is  only  valid  until  the  year  2050  due  to  truncation  (see  Vallado  pg  265) . 


7.  note  all  of  these  values  must  be  in  degrees,  not  radians!  !  ! 


JD  =  juldatet(t);  %  JD  is  in  ref  to  UTC  (let  JD_UTC  =  JD_UT1) 

7.  find  the  number  of  Julian  centuries 

T_UT1  =  (JD  -  2451545. 0)/36525;  7.  Vallado  pg  188 

lambda_M_sundeg  =  280.460  +  36000 . 770*T_UT1 ;  7.  deg 

T_TDB  =  T_UT1 ; 

M_sundeg  =  357.5277233  +  35999 . 05034*T_TDB ;  7.  deg 

lambda_ecliptic  =  lambda_M_sundeg  +  1 . 914666471*sind(M_sundeg)  + 
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0 . 019994643*sind(2*M_sundeg) ;  %  deg 

r_sim  =  1.000140612  -  0 . 016708617*cosd(M_sundeg) 

0 . 000139589*cosd(2*M_sundeg) ;  %  AU 
epsilon  =  23.439291  -  0 . 0130042*T_TDB ;  7.  deg 

7,  the  vector  referred  to  the  epoch  of  date  is  (in  AU) 
r_ES_vecAU  =  [r_sun*cosd(lambda_ecliptic) ; . . . 

r_sun*cosd(epsilon) *sind(lambda_ecliptic) ; . . . 
r_sun*sind (epsilon) *sind(lambda_ecliptic)]  ; 

r_ES_vec  =  149598000 . *r_ES_vecAU;  7o  Sun-Earth  vector  in  km 

7.  also  calculate  the  right  ascension  and  declination  of  the  sun  (can  check 
7.  with  Astronomical  Almanac) 

alpha_sun  =  atand (cosd (epsilon) *tand (lambda_ecliptic) ) ;  delta_sun 
=  asind(sind(epsilon)*sind(lambda_ecliptic)) ; 

return 


C-4  EQUcode.m 

7.  This  m-file  is  called  by  the  executable,  EQUcode_wrapper .m 


7.  Define  initial  orbital  parameters  for  the  DSCS  II  F-13  geostationary 
7.  satellite  based  on  a  two  line  element  set  at  specific  UTC  (epoch  at  yr  05 
7.  day  202  fraction  .40912884) 


7.  Line  2  of  TLE : 

BC_o  =  -.00000117;  %  ballistic  coefficient 

epoch_yr  =  2005;  7»  year  in  4-digit  format 

epoch_day  =  202;  epoch_frac  =  .40912884; 

JD_o  =  juldate(epoch_yr ,epoch_day,epoch_frac) ; 


7»  inclination  (deg) 
%  inclination  (rad) 


7.  Line  3  of  TLE: 
i_o_deg  =  14.4997; 

i_o  =  i_o_deg*2*pi()/360; 

RAAN_o_deg  =  18.8800; 

7.  right  ascension  of  the  ascending  node  (deg) 

RAAN_o  =  RAAN_o_deg*2*pi () /360 ; 

7.  right  ascension  of  the  ascending  node  (rad) 

7.  eccentricity  (dimensionless) 
l  argument  of  perigee  (deg) 
omega_o  =  omega_o_deg*2*pi () /360 ;  7«  argument  of  perigee  (rad) 

M_o_deg  =  38.4719; 

M_o  =  M_o_deg*2*pi()/360; 
n_o_rev  =  .95598989; 

n_o  =  n_o_rev* (2*pi () /86400) ; 
revnum.o  =  54242 ; 


e_o  =  0.0005671; 
omega_o_deg  =  321.5388; 


7»  mean  anomaly  (deg) 

7.  mean  anomaly  (rad) 

7.  mean  motion  (rev/day) 

7.  mean  motion  (rad/sec) 

7»  revolution  number  at  epoch  (rev) 
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%  Calculate  initial  parameters : 


a_o  =  (mu_earth/n_o~2) ~ (1/3) ; 

ra_o  =  a_o*(l  +  e_o) ; 
rp_o  =  a_o*(l  -  e_o) ; 

ha_o  =  ra_o  -  R_earth; 
hp_o  =  rp_o  -  R_earth; 

T_o  =  2*pi () *sqrt (a_o~3/mu_earth) ; 

°/0  approximate  period  of  orbit  (sec)  ,  Hale  pg  62  (unperturbed) 

T_o_hr  =  T_o/ (60*60)  ; 

%  approximate  period  of  orbit  (hr)  (unperturbed) 

Energy_o  =  -mu_earth/ (2*a_o) ; 

°/0  total  specific  mechanical  energy  (km~2/sec~2)  ,  Hale  pg  18 
t_o  =  (T_o/(2*pi () ) ) *M_o ; 

%  time  elapsed  from  perigee  passage  (sec) ,  Hale  pg  63 

%  Eccentric  anomaly  evaluation  using  Newton-Raphson  iteration 
E_o  =  CalcEA(M_o,e_o, le-10) ; 

nu_o  =  2*atan(sqrt ( (l+e_o) /(l-e_o) ) *tan(E_o/2) ) ;  %  true  anomaly  (rad),  Hale 
r_o  =  a_o*(l  -  e_o*cos (E_o) ) ;  70  instantaneous  radius  of  s/c  (km),  Hale 

h_o  =  r_o  -  R_earth;  7,  instantaneous  altitude  of  s/c  (km) 

rho  =  rho_s*exp (-beta*h_o)  ;  °L  density  at  h_o  (kg/knT3) 


7.  semimajor  axis  (km),  Hale  pg  66 

7.  radius  of  apoapsis  (km)  ,  Hale 
7.  radius  of  periapsis  (km)  ,  Hale 


7.  height  of  apoapsis  (km)  ,  Hale 
7.  height  of  periapsis  (km)  ,  Hale 


7,  Calculate  the  geocentric  position  vector  of  the  Sun: 


r_ES_vec  =  earthsun(epoch_yr,epoch_day,epoch_frac) ;  %  returns  in  units  of  km 


7.  Calculate  the  initial  modified  equinoctial  (MEE)  element  set: 


p_o  =  a_o*(l  -  e_o~2)  ;  %  semilatus  rectum  (km) 

f_o  =  e_o*cos (omega_o  +  RAAN_o) ;  g_o  =  e_o*sin(omega_o  +  RAAN_o) ; 
h_o  =  tan(i_o/2)*cos(RAAN_o) ;  k_o  =  tan(i_o/2)*sin(RAAN_o) ; 

L_o  =  RAAN_o  +  omega_o  +  nu_o;  7«  true  longitude 


7.  Call  function  mequtorv  to  convert  the  initial  modified  MEE 
7.  element  set  to  inital  r  and  v  vectors 

x_o_vec  =  mequtorv(p_o,  f_o,  g_o,  h_o,  k_o,  L_o) ; 
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%  of  the  form  [xo;  yo;  zo;  xdoto;  ydoto;  zdoto]  with  r  components  in  km 
%  and  v  components  in  km/sec 

r_o  =  [x_o_vec(l);  x_o_vec(2);  x_o_vec(3)]  ;  v_o  =  [x_o_vec(4); 
x_o_vec(5);  x_o_vec(6)]  ; 


%  Call  function  propxequ  to  propagate  the  state  vector  x  from  x_o  using 
1  the  relationship  xdot  =  f(x,t) 

aa  =  1 ;  equo_vec  =  [p_o;  f_o;  g_o;  h_o;  k_o;  L_o]  ; 

ti  =  0;  options=odeset ( 'MaxStep 1 ,MaxStep) ;  [t,x]  = 
ode45 (Opropxequ,  [ti  n*T_o] , equo_vec , options) ; 

%  [t,x]  =  ode45 (Opropxequ, tspan,equo_vec, options) ;  7.  uncomment  to  output  a 
7,  vector  with  specific  time  intervals;  use  with  STK  results  to  plot  error 

poft  =  x(:,l);  foft  =  x(:,2);  goft  =  x(:,3);  hoft  =  x(:,4);  koft  = 
x(:,5);  Loft  =  x(:,6); 


7.  compute  r(t)  and  v(t)  based  on  the  MEE  elements 


for  i=l : length (poft) 

alpha  =  sqrt (hoft (i) ~2  -  koft(i)~2); 
s  =  sqrt(l  +  hoft(i)~2  +  koft(i)"2); 
w  =  1  +  f oft (i) *cos (Loft  (i) )  +  goft(i)*sin(Loft(i)) ; 
r  =  poft(i)/w; 

rl  =  (r/s~2) * (cos (Loft (i) )  +  alpha~2*cos (Loft (i) )  +... 

2*hoft (i) *koft (i) *s in (Loft (i) ) ) ; 
r2  =  (r/s~2) * (sin(Loft (i) )  -  alpha~2*sin(Loft (i) )  +... 

2*hoft (i) *koft (i) *cos (Loft  (i) ) )  ; 
r3  =  (2*r/s~2) * (hoft (i) *sin(Loft (i) )  -  koft (i) *cos (Loft (i) )) ; 
r  =  [rl  r2  r3] ; 

vl  =  (-l/s~2) *sqrt (mu_earth/poft (i) ) * (sin(Loft (i) )  +  alpha~2*sin(Loft (i) ) . 
-  2*hoft (i) *koft (i) *cos (Loft (i) )  +  goft(i)  -  2*f oft (i) *hoft (i) *koft (i) 
+  alpha~2*goft (i) ) ; 

v2  =  (-l/s~2) *sqrt (mu_earth/poft (i) )* (-cos (Loft (i) )  +  alpha~2*cos (Loft (i) ) 
+  2*hoft (i) *koft (i) *sin(Loft (i) )  -  foft(i)  +  2*goft (i) *hoft (i) *koft (i) 
+  alpha~2*f oft (i) )  ; 

v3  =  (2/s~2)*sqrt(mu_earth/poft(i))*(hoft(i)*cos(Loft(i))  +... 

koft(i)*sin(Loft(i))  +  f oft (i) *hoft (i)  +  goft (i) *koft (i) ) ; 
v  =  [vl  v2  v3] ; 

rvec  =  [rvec ;  r]  ; 
wee  =  [wee ;  v]  ; 

end 
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°/0  convert  MEE  set  back  to  Keplerian  element  set : 


for  ii  =  l:length(x) 

a  =  poft(ii)/(l  -  foft(ii)~2  -  goft  (ii)  ~2)  ;  7o  semimajor  axis  (km) 
a_vec  =  [a_vec  a] ; 

e  =  sqrt ( (f oft (ii) ) ~2  +  (goft (ii) ) ~2) ;  %  eccentricity  (dimensionless) 
e_vec  =  [e_vec  e] ; 

i  =  atan2 (2*sqrt ( (hoft (ii) ) ~2  +  (koft (ii) ) ~2) , 1-hoft (ii) ~2-koft (ii) ~2)  ; 

7.  inclination  (rad) 
i_vec  =  [i_vec  i] ;  %  rad 

w  =  atan2 (goft (ii) *hoft (ii) -f oft (ii) *koft (ii) ,foft(ii)*hoft(ii)+. . . 
goft (ii) *koft (ii) ) ; 

7.  argument  of  periapsis  (rad) 
w  =  rad2deg(w)  ;  °/0  deg 
w  =  zero22pi(w);  7.  deg 
w_vec  =  [w_vec  w]  ;  7o  deg 

RAAN  =  atan2 (koft (ii) ,hoft(ii)) ; 

7.  right  ascension  of  the  ascending  node  (rad) 

RAAN  =  rad2deg(RAAN) ;  /  deg 
RAAN_vec  =  [RAAN_vec  RAAN]  ;  °L  deg 

mi  =  Loft(ii)  -  (atan2 (koft (ii) , hoft (ii) )  +... 

atan2 (goft (ii) *hoft (ii) -f oft (ii) *koft (ii) ,foft(ii)*hoft(ii)+. . . 
goft  (ii)  *koft  (ii)  )  )  ;  7o  true  anomaly  (rad) 
nu.  =  rad2deg(nu)  ;  7o  deg 
mi  =  zero22pi  (nu)  ;  7o  deg 
nu_vec  =  [nu_vec  nu]  ;  7o  deg 

u  =  atan2 (hoft (ii) *sin(Loft (ii) ) -koft (ii) *cos (Loft (ii) ),..  . 
hof t (ii) *cos (Loft (ii) )+koft (ii) *s in (Loft (ii) ) ) ; 
l  argument  of  latitude 
u_vec  =  [u_vec  u] ; 

end 

i_vec_deg  =  i_vec  .  *180/pi  ()  ;  °L  deg 

w_vec_deg  =  w_vec; 


7.  Calculate  misc.  parameters  of  interest 


rp_vec  =  a_vec.*(l  -  e_vec)  ;  7«  perigee  altitude  (km) 

OrbitRaise  =  a_vec(length(a_vec))-a_vec(l) ;  a_end  = 
a_vec (length(a_vec) ) ; 
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elapsedtime_days  =  n*T_o/60/60/24;  elapsedtime_years  = 
n*T_o/ 60/ 60/24/365; 


omega_i  =  w_vec_deg(l) ;  omega_f  =  w_vec_deg(end) ; 

i_i  =  i_vec_deg(l) ;  i_f  =  i_vec_deg(end) ; 

RAAN_i  =  RAAN_vec (1) ;  RAAN_f  =  RAAN_vec(end) ; 

xx  =  polyf  it  (t  ’  ,  a_vec ,  1)  ;  polx  =  xx(l)  *t+xx(2)  ;  a_ch.an.ge  = 
polx(length(polx) )-polx(l) ; 

%  plot : 


%  uncomment  figs  1-3  for  baseline  case  plots 

7.  figured) 

7.  subplot (2,1,1); 

%  plot(t,rvec(: ,1) , *b: ’) 

7.  hold  on 

7.  plot  (t  ,rvec(  :  ,2)  ,  ’g — ’) 

7,  hold  on 

7.  plot  (t  ,rvec(  :  ,3)  ,  ’r’  ) 

7.  title ([ ’  Position  vector  versus  time  over  ’ ,num2str (n) , ’  orbital  periods1]) 
7.  legend  ( 1  rx 1  ,  1  ry 1 ,  ’  rz 1 ) 

7.  xlabelC’time  (sec)’) 

7.  ylabel( ’position  (km)’) 

7.  subplot  (2d  ,2) ; 

7,  plot (t, vvec(:  d)  db:  ’) 

7.  hold  on 

7.  plot  (t,  wee  (  :  ,2)  dg — ’  ) 

7.  hold  on 

7.  plot  (t ,  vvec(  :  ,3)  dr1 ) 

7,  title( [’Velocity  vector  versus  time  over  ’ ,num2str (n) d  orbital  periods’]) 
7.  legend  (  ’  vx  ’  d  vy  ’  d  vz  ’ ) 

7.  xlabelC’time  (sec)’) 

7.  ylabel( ’velocity  (km/s)’) 

7.  print  (1,  ’-depsc’,  ’randvl’) 

7. 

7.  figure (2) 

7.  plot(rvec(:  d)  ,rvec(:  ,2)) 

7,  axis  square 

7.  xlim(  [-50000  50000]) 

7.  ylim(  [-50000  50000]) 

7.  hold  on 

7,  ttt=linspace(0,2*pi)  ! 

7.  f ill(6378*cos(ttt)  ,6378*sin(ttt)  ,[.5  .8  .3]) 

7.  axis  square 
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%  hold  on 

7.  plot(rvec(l, 1) ,rvec(l,2) , ’b-.o’ , ’LineWidth’ ,2, ’MarkerEdgeColor 1 , ’r’ , . . . 

7.  ’MarkerFaceColor ’ ,’r’, ’MarkerSize’ ,5) 

°/  text(rvec(l , 1) ,rvec(l,2) , 1 \lef tarrowStarting  value’) 

%  axis  square 

%  title([’2-D  orbit  trace  of  satellite  over  ’ ,num2str (n) , ’  orbital  periods’]) 
'/,  xlabel(’rx  (km)’) 

'/,  ylabel(’ry  (km)’) 

1  hold  off 

%  print (2,  ’-depsc’,  ’2-Dtracel’) 

7. 

7.  figure (3) 

7.  plot3  (rvec(  :  ,  1)  ,rvec(  :  ,  2)  ,rvec(  :  ,3)  ) 

7,  title([’3-D  orbit  trace  of  satellite  over  ’ ,num2str (n) , ’  orbital  periods’]) 
7.  grid  on 
7.  axis  tight 
7.  xlim(  [-40000  40000]) 

7.  ylim(  [-40000  40000]) 

7.  zlim(  [-40000  40000]) 

7.  hold  on 

7.  colormap  winter 

7.  [X,Y,Z]  =  sphere (20); 

7.  surf  (6378*X,6378*Y,6378*Z) 

7.  xlabel(’rx  (km)’) 

7.  ylabel(’ry  (km)’) 

7.  zlabel(’rz  (km)’) 

7.  print (3,  ’-depsc’,  ’3-Dtracel’) 

switch  nn 

case  1  7.  1  orbit 

meanval  =  mean(a_vec) ; 
if  meanval>a_vec(l) 

achange  =  meanval  -  a_vec(l); 

else 

achange  =  a_vec(l)  -  meanval; 

end 

slope  =  achange/t(length(t)) ; 
figure (4) ; 

subplot (3,1,1);  plot (t , i_vec_deg, ’ r- ’ ) 
title([’\it  i  \rm  versus  time  over  ’ ,num2str(n) , . . . 

’  orbital  periods ’],’ fontsize ’, 12) 
set (gca, ’fontsize ’ , 12) 
xlabel( ’time  (sec) ’ , ’fontsize’ , 12) 
ylabel(’\it  i  \rm  (deg) ’,’ fontsize ’, 12) 
subplot (3,1,2);  plot (t , a_vec , ’b- ’ ) 
set (gca, ’fontsize ’ , 12) 

title([’\it  a  \rm  versus  time  over  ’ ,num2str(n) , . . . 

’  orbital  periods;  \Delta\it  a  \rm=’ ,num2str (achange) , ’  km’],... 
’fontsize’ , 12) 

xlabel( ’time  (sec) ’ , ’fontsize’ , 12) 
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ylabel (’\it  a  \rm  (km) 1 , ’ f ontsize 1 , 12) 
subplot (3,1,3);  plot (t , e_vec , ’ g- ’ ) 
set (gca, ’  font size ’  ,  12) 

title([’\it  e  \rm  versus  time  over  ’ ,num2str (n) , ’  orbital  periods1] 
’fontsize’ , 12) 

7,  ylim(  [0  0.002]) 

xlabel( ’time  (sec) 1 , ’fontsize’ , 12) 
ylabel(’\it  e ’,’ fontsize ’, 12) 

end 

switch  nn 

case  {2,3,4} 
figure (4) ; 

subplot (3,1,1);  plot (t , i_vec_deg, ’ r- ’ ) 
title([’\it  i  \rm  versus  time  over  ’ ,num2str (n) , . . . 

1  orbital  periods ’],’ fontsize ’, 12) 
set (gca, ’ font size ’ , 12) 
xlabel( ’time  (sec) ’ , ’ font size ’ , 12) 
ylabel(’\it  i  \rm  (deg) ’,’ f ontsize ’, 12) 
subplot (3,1,2);  plot (t , a_vec , ’b- ’ ) 
set (gca, ’ font size ’ , 12) 

title([’\it  a  \rm  versus  time  over  ’ ,num2str (n) , . . . 

’  orbital  periods;  \Delta\it  a  \rm=’ ,num2str (a_change) , ’  km’],, 
’fontsize’ , 12) 

xlabel( ’time  (sec) ’ , ’fontsize’ , 12) 
ylabel(’\it  a  \rm  (km) ’,’ fontsize ’, 12) 
subplot (3,1,3);  plot (t , e_vec , ’ g- ’ ) 
set (gca, ’fontsize ’ , 12) 

title([’\it  e  \rm  versus  time  over  ’ ,num2str (n) , ’  orbital  periods’] 
’fontsize’ , 12) 

*/.  ylim(  [0  0.002]) 

xlabel ( ’time  (sec) ’ , ’fontsize ’ , 12) 
ylabel(’\it  e ’,’ fontsize ’, 12) 

end 

figure (5) ; 

subplot (3, 1 , 1)  ;  plot (t , w_vec_deg, ’g- ’ ) ;  title ([’ \omega  versus  time 
over  ’ ,num2str (n) , ’  orbital  periods’] , ’fontsize’ , 12) ; 
set (gca, ’ fontsize ’, 12) ;  xlabel(’time  (sec) ’,’ fontsize ’, 12) ; 
ylabel ( ’ \omega  (deg) ’ , ’fontsize ’ , 12) 

subplot (3, 1 , 2) ;  plot (t ,RAAN_vec , ’b- ’ ) ;  title ( [’\0mega  versus  time 
over  ’ ,num2str (n) , ’  orbital  periods’] , ’fontsize’ , 12) ; 
set (gca, ’ fontsize ’, 12) ;  xlabel(’time  (sec) ’,’ fontsize ’, 12) ; 
ylabel ( ’ \0mega  (deg) ’ , ’fontsize ’ , 12) 

subplot (3, 1 , 3) ;  plot(t,nu_vec, ’k-’) ;  title ([’\nu  versus  time  over 
’ ,num2str (n) , ’  orbital  periods ’],’ fontsize ’, 12) ; 
set (gca, ’ fontsize ’, 12) ;  xlabel(’time  (sec) ’,’ fontsize ’, 12) ; 
ylabel ( ’ \nu  (deg) ’ , ’fontsize’ ,12) 
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C.5  EQUcode-wrapper.m 

V.  This  script  file  is  the  executable  portion  of  the  thesis  code.  It  calls 
%  EQUcode.m,  and  includes  switch-case  functionality  to  enable  "one-push  of 
%  the  button"results .  Recommend  running  through  Linux,  as  the  simulation 
%  takes  approximately  12  hours  to  push  through  all  data  sets.  The 
%  simulation  saves  off  pertinanent  workspace  variables  and  all  plots  (eps)  . 

%  Two-body  orbit  determination  code  in  modified  equinoctial  elements 
7,  Capt  Jody  Paris  (GA-06M) 

clc 

format  long 
syms  E 

a_vec  =  []  ;  e_vec  =  []  ;  i_vec  =  []  ;  w_vec  =  []  ;  RAAN_vec  =  []  ; 
nu_vec  =  []  ;  u_vec  =  []  ;  rvec  =  []  ;  wee  =  []  ; 


%  Define  constants  and  global  variables 


global  n  g_s  G  M_earth  R_earth  mu_earth  mu_moon  rho  J2  SF  c  m 
epoch_yr . . . 

epoch_day  epoch_frac  A_plate  A_cyl  Earth_aphelion_table  c_R_plate . . . 
c_R_cyl  aa  s_vec  cont  moon  countt 


g_s  =  9 . 796432222e-3 ; 

%  mean  gravitational  acceleration  of  Earth  at  surface  (km/s~2) 


G  =  6 . 67259e-ll ; 
M_earth  =  5.9T3Te24; 
R_earth  =  6378; 
mu_earth  =  3.986004e5; 
beta  =  0.14; 
rho_s  =  1.225e9; 

J2  =  0.0010826269; 

SF  =  1353; 
c  =  299792458/1000; 
m_moon  =  7.3483e22; 
mu_moon  =  4902.799; 
Earth_aphelion_table  = 


7.  universal  gravitational  constant  (Nm~2/kg~2) 
l  Earth’s  mass  (kg) 

7.  Earth’s  radius  (km) 

7.  Earth’s  gravitational  parameter  (km~3/s~2) 

7.  km~-l 

7.  density  at  surface  of  Earth,  kg/km~3 
7.  Earth’s  oblateness  coefficient 
7.  solar  flux  value  (kg/sec~2) 

7.  speed  of  light  (km/sec) 

7.  mass  of  moon  (kg) 

7,  gravitational  parameter  of  moon  (km~3/s~2) 
[2005  186 ; 2006  184;2007  188;2008  186;2009 


185; .. . 

2010  187 ; 2011  185;2012  187;2013  186;2014  185;2015  187;2016  186;... 
2017  184; 2018  187;2019  185;2020  186]; 

7.  in  [year  date]  format,  from  US  Naval  Observatory 


m  =  611; 

1  =  .004; 
d  =  .003; 

A_plate  =  .00000335; 


7,  satellite  mass,  kg 

7.  satellite  length  (rotating  platform)  ,  km 
7.  satellite  diameter,  km 

7.  area  of  flat  plate  (parabolic  antennae)  ,  km~2 
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*/.  area  of  spinning  body  (cylindrical)  ,  km~2 


A_cyl  =  .00000502; 
c_R_cyl  =  [1/3;  1/3;  1/3]  ; 

7.  reflectivity  coefficients  [absorption;  diffuse;  specular]  for 
V.  cylinder 


%  run  baseline  case  (no  controller,  despun  area  =  0) 

cont  =  1;  controller  =  ’b’; 
c_R_plate  =  [1/3;  1/3;  1/3]; 

%  reflectivity  coefficients  [absorption;  diffuse;  specular] 
for  nn=l : 3 

switch  nn 
case  1 

n  =  1 ;  7.  number  of  orbits  to  propagate  over 

MaxStep  =  le2; 
case  2 

n  =  100;  7.  number  of  orbits  to  propagate  over 

MaxStep  =  le3; 
case  3 

n  =  1000;  7,  number  of  orbits  to  propagate  over 

MaxStep  =  le4; 
case  4 

n  =  10000; 

MaxStep  =  le5; 

end 

multf actor  =  1; 

filename  =  [’b_J  num2str(n)]; 


y 

i 

7. 

7. 


7. 

7. 

7. 

7. 

7. 

7. 

7. 


EQUcode  7«  call  EQUcode.m 
7.  help  save 

filenamet  =  [filename  ’t’]; 
filenamei  =  [filename  ’i’]; 
filenamea  =  [filename  ’a’]; 
filenamee  =  [filename  ’e’]; 


filenamew  =  [filename  ’w’]; 

7,  uncomment 

to 

save 

filenameRAAN  = 

[filename  ’ RAAN’]; 

7,  uncomment 

to 

save 

filenamenu  =  [filename  ’nu’] ; 

7.  uncomment 

to 

save 

save (f ilenamet , ’t 

’) 

save (f ilenamei , ’ i 

_vec_deg’ ) 

save (f ilenamea, ’a. 

_vec ’ ) 

save (f ilenamee , ’e. 

_vec ’ ) 

save (filenamew , 

1  w_vec_deg’ ) 

7.  uncomment 

to 

save 

save (filenameRAAN, ’RAAN_vec’) 

7.  uncomment 

to 

save 

save (filenamenu 

,  ’nu_vec ’ ) 

7,  uncomment 

to 

save 

save( ’poftlOOO1 

,  1 pof t ’ ) 

7.  uncomment 

to 

save 

save ( ’ gof tlOOO ’ 

,  ’goft ’ ) 

7,  uncomment 

to 

save 

save ( 1 f oftlOOO 1 

,  ’foft’) 

7,  uncomment 

to 

save 

save( ’LoftlOOO1 

,  1  Loft 1 ) 

7,  uncomment 

to 

save 
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%  save ( ’koftlOOO 1 , ’koft ’ ) 
°/0  save  (  ’hoft  1000 1  ,  ’hoft  ’  ) 


%  uncomment  to  save 
%  uncomment  to  save 


V.  help  print 

filenamef4  =  [filename  ’iae’]; 
filenamef5  =  [filename  ’wRAANnu’]; 
print(4,  ’-depsc’,  filenamef4) 
print(5,  ’-depsc’,  filenamef5) 


%  prints  figure  4  to  .eps  file 
V.  prints  figure  5  to  .eps  file 


close  all; 
clear  filename 

a_vec  =  []  ;  e_vec  =  []  ;  i_vec  =  []  ;  w_vec  =  [];... 

RAAN_vec  =  []  ;  nu_vec  =  []  ;  u_vec  =  []  ;  rvec  =  [];... 
wee  =  []  ; 


end 


%  run  other  cases 

for  cont=2:3  %  contoller  is  indexed  in  srpaccel.m 
switch  cont 
case  1 


controller  = 

’b’; 

1  baseline  case 

case  2 

controller  = 
cc  =  ’simple1 

’s’  ; 

> 

*/,  simple  controller 

case  3 

controller  = 

’c’  ; 

*/,  complex  controller 

cc  =  1  complex’ ; 

end 

for  CRP=1:5 

switch  CRP 
case  1 

c_R_plate  =  [1/3;  1/3;  1/3] ; 

7.  reflectivity  coefficients 

7  [absorption;  diffuse;  specular]  for  antennae 
coeff  =  [controller  ’-thirds’]; 
case  2 

c_R_plate  =  [1;  0;  0];  °/»  all  absorption 

coeff  =  [controller  ’_a’]; 
case  3 

c_R_plate  =  [0;  0;  1];  °/»  all  specular  reflection 

coeff  =  [controller  ’_s’]; 
case  4 

c_R_plate  =  [0;  1;  0];  l  all  diffuse  reflection 
coeff  =  [controller  ’_d’] ; 
case  5 

c_R_plate  =  [0;  .5;  .5];  °i  specular  &  diffuse  reflection 
coeff  =  [controller  ’_sdJ]; 

end 
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for  nn=l : 3 

7.  Define  varying  parameters  for  code 

switch  nn 
case  1 

n  =  1 ;  7,  number  of  orbits  to  propagate  over 

MaxStep  =  le2; 
orbitt  =  [coeff  ’_1’]; 
case  2 

n  =  100;  7.  number  of  orbits  to  propagate  over 

MaxStep  =  le3; 
orbitt  =  [coeff  *  _100  * ] ; 
case  3 

n  =  1000;  7.  number  of  orbits  to  propagate  over 

MaxStep  =  le4; 

orbitt  =  [coeff  3  _ 1000 J ] ; 

end 

for  mf =1 : 3 

switch  mf 
case  1 

multf actor  =  1; 

A_plate  =  . 00000335*multf actor ; 
filename  =  [orbitt  3_13]; 
case  2 

multf actor  =  5; 

A_plate  =  . 00000335*multf actor ; 
filename  =  [orbitt  3 _5 3 ] ; 
case  3 

multf actor  =  10; 

A_plate  =  . 00000335*multf actor ; 
filename  =  [orbitt  3 _ 10 3 ]  ; 

end 

EQUcode 


7.  help  save 


f ilenamet  =  [filename 

]; 

filenamei  =  [filename 

;  i 

]; 

filenamea  =  [filename 

’  a 

]; 

filenamee  =  [filename 

J  e 

]; 

filenamew  =  [filename 

;  w 

]; 

filenameRAAN  =  [filename  1 RAAN1]; 
filenamenu  =  [filename  ’nu’] ; 

save (f ilenamet , 1 t ’ ) 
save (f ilenamei , 1 i_vec_deg’ ) 
save (f ilenamea, ,a_vec)) 
save (f ilenamee , ,e_vec)) 
save (f ilenamew , 1 w_vec_deg ’ ) 
save (filenameRAAN, ,RAAN_vec3) 
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save (f ilenamenu, ’nu_vec’) 


7,  help  print 

filenamef4  =  [filename  ’iae’]; 
filenamef5  =  [filename  ’wRAANnu’]; 
print(4,  ’-depsc’,  filenamef4) 
print(5,  ’-depsc’,  filenamef5) 


if  cont~=l 

bt  =  [’b_’  num2str(n)  ’t’]; 
tb=load(bt) ; 

°/„  baseline  time  vector  for  current  case 
ba  =  [’b_’  num2str(n)  ’a’]; 
ab=load(ba) ; 

7„  baseline  semimajor  axis  vector  for  current  case 

x  =  polyf it (t ’ , a_vec , 1)  ; 
polx  =  x(l)*t+x(2); 

a_change  =  polx (length (polx) )-polx(l)  ; 
figure (6) 

plot(tb.t,ab.a_vec, ’b — ’) 

title([’\it  a  \rm  versus  time  over  1 ,num2str (n) , . . . 

’  orbital  periods;  \Delta\it  a  \rm=I , . . . 
num2str (a_change) ,  ’  km’] , ’fontsize’ , 12) 
set (gca, ’fontsize ’ , 12) 
xlabel ( ’time  (sec) ’ , ’fontsize ’ , 12) 
ylabel(’\it  a  \rm  (km) ’,’ fontsize ’, 12) 
hold  on 

plot (t , a_vec , ’ r ’ ) 

legend ( ’baseline ’ ,num2str (cc) , ’ Location’ , ’Best ’ ) 
hold  off 

filenamef6  =  [’mix_’  filename]; 
print(6,  ’-depsc’,  filenamef6) 


end 

close  all 
clear  filename 

a_vec  =  []  ;  e_vec  =  []  ;  i_vec  =  []  ;  w_vec  =  [];... 

RAAN_vec  =  []  ;  nu_vec  =  []  ;  u_vec  =  []  ;  rvec  =  [];... 
wee  =  []  ; 


end 


end 


end 
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end 


C.6  juldate.m 

function  JD  =  juldate(epoch_yr,epoch_day,epoch_frac) 

%  This  function  inputs  the  epoch  information  given  in  the  current  TLE  set 

%  and  outputs  the  appropriate  Julian  day  based  on  the  algorithm  given  on 

°/0  page  186  of  Vallado.  Here,  JD  is  based  on  UTC  but  can  be  set  equal  to  JD_UT1. 

%  calculate  the  following  based  on  epoch_frac  and  epoch_day: 
if  epoch_yr/4  ==  f ix(epoch_yr/4) 

if  epoch_yr/100  ==  fix(epoch_yr/100) 

if  epoch_yr/400  ==  f ix(epoch_yr/400) 


j  =  i; 

1  leap  year 

else 

j  =  0; 

l  not  a  leap  year 

end 

j  =  i; 


end 

else 

j  =  0; 

end 


if  epoch_day  <=  31 

epoch_mo  =  1; 
epoch_dofm  =  epoch_day; 
elseif  epoch_day  <=  59+j 
epoch_mo  =  2; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  90+j 
epoch_mo  =  3; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  120+j 
epoch_mo  =  4; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  151+ j 
epoch_mo  =  5; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  181+ j 
epoch_mo  =  6; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  212+j 
epoch_mo  =  7 ; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  243+j 
epoch_mo  =  8; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  273+j 
epoch_mo  =  9; 


-  31+ j ; 

-  59+j ; 

-  90+j ; 


-  120+j; 


-  151+ j  ; 


-  181+ j ; 


-  212+j; 
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epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  304+j 
epoch_mo  =  10; 
epoch_dofm  =  epoch_day 
elseif  epoch_day  <=  334+j 
epoch_mo  =  11; 
epoch_dofm  =  epoch_day 

else 

epoch_mo  =  12; 
epoch_dofm  =  epoch_day 

end 


243+ j  ; 


273+j ; 


304+j ; 


334+j ; 


epoch_hrlong  =  epoch_frac*24;  7.  chop  off  remainder  and  use  for  min 
epoch_hr  =  f ix(epoch_hrlong) ; 
epoch_minlong  =  (epoch_hrlong  -  epoch_hr) *60 ; 

7,  chop  off  remainder  and  use  for  sec 
epoch_min  =  f ix(epoch_minlong) ; 
epoch_sec  =  (epoch_minlong  -  epoch_min) *60 ; 

JD  =  367*epoch_yr-real((7*epoch_yr+real((epoch_mo+9)/12))/4)+. . . 
real (275*epoch_mo/9) +epoch_dofm+1721013 . 5+ ( ( ( (epoch_sec/ 60) + . . . 
epoch_min) / 60) +epoch_hr ) /24 ; 

return 


C.  1  juldatet.m 

function  JD  =  juldatet(t) 

7.  This  function  inputs  the  epoch  information  given  in  the  current  TLE  set 

7,  and  outputs  the  appropriate  Julian  day  based  on  the  algorithm  given  on 

7.  page  186  of  Vallado.  Here,  JD  is  based  on  UTC  but  can  be  set  equal  to  JD_UT1. 


epoch_yr_o  =  2005;  7,  year  in  4-digit  format 
epoch_day_o  =  202;  epoch_frac_o  =  .40912884; 

current_epoch_yrlong  =  t*3. 1688764640818e-8;  7«  convert  t  in  sec  to  years 
current_epoch_yr  =  f ix(current_epoch_yrlong) ; 

current_remainder_yr  =  current_epoch_yrlong  -  f ix(current_epoch_yrlong) ; 
current_epoch_daylong  =  current_remainder_yr*365 . 2442 ; 

7.  remainder  of  years  to  days 
current_epoch_day  =  f ix(current_epoch_daylong) ; 

current_epoch_frac  =  current_epoch_daylong  -  f ix(current_epoch_daylong) ; 
7.  fraction  of  a  day  leftover 

epoch_yr  =  epoch_yr_o  +  current_epoch_yr ;  epoch_day  =  epoch_day_o 
+  current_epoch_day ;  epoch_frac  =  epoch_frac_o  + 
current_epoch_f rac ; 

7.  calculate  the  following  based  on  epoch_frac  and  epoch_day: 
if  epoch_yr/4  ==  f ix(epoch_yr/4) 
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if  epoch_yr/100  ==  fix(epoch_yr/100) 

if  epoch_yr/400  ==  f ix(epoch_yr/400) 


j  =  i; 

1  leap  year 

else 

j  =  0; 

%  not  a  leap  year 

end 

j  =  i; 


end 

else 

j  =  0; 

end 

if  epoch_day  <=  31 

epoch_mo  =  1; 
epoch_dofm  =  epoch_day; 
elseif  epoch_day  <=  59+j 
epoch_mo  =  2; 

epoch_dofm  =  epoch_day  -  31+j ; 
elseif  epoch_day  <=  90+j 
epoch_mo  =  3; 

epoch_dofm  =  epoch_day  -  59+j ; 
elseif  epoch_day  <=  120+j 
epoch_mo  =  4; 

epoch_dofm  =  epoch_day  -  90+j ; 
elseif  epoch_day  <=  151+ j 
epoch_mo  =  5; 

epoch_dofm  =  epoch_day  -  120+j ; 
elseif  epoch_day  <=  181+ j 
epoch_mo  =  6; 

epoch_dofm  =  epoch_day  -  151+ j ; 
elseif  epoch_day  <=  212+j 
epoch_mo  =  7 ; 

epoch_dofm  =  epoch_day  -  181+ j ; 
elseif  epoch_day  <=  243+j 
epoch_mo  =  8; 

epoch_dofm  =  epoch_day  -  212+j ; 
elseif  epoch_day  <=  273+j 
epoch_mo  =  9; 

epoch_dofm  =  epoch_day  -  243+j ; 
elseif  epoch_day  <=  304+j 
epoch_mo  =  10; 

epoch_dofm  =  epoch_day  -  273+j ; 
elseif  epoch_day  <=  334+j 
epoch_mo  =  11; 

epoch_dofm  =  epoch_day  -  304+j ; 

else 

epoch_mo  =  12; 

epoch_dofm  =  epoch_day  -  334+j ; 

end 
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epoch_hrlong  =  epoch_frac*24;  7.  chop  off  remainder  and  use  for  min 
epoch_hr  =  f ix(epoch_hrlong) ; 
epoch_minlong  =  (epoch_hrlong  -  epoch_hr) *60 ; 

7,  chop  off  remainder  and  use  for  sec 
epoch_min  =  f ix(epoch_minlong) ; 
epoch_sec  =  (epoch_minlong  -  epoch_min) *60 ; 

JD  =  367*epoch_yr-real((7*epoch_yr+real((epoch_mo+9)/12))/4)+. . . 
real (275*epoch_mo/9) +epoch_dofm+1721013 . 5+ ( ( ( (epoch_sec/ 60) + . . . 
epoch_min) / 60) +epoch_hr ) /24 ; 

return 


C.8  mequtorv.m 

function  x_o_vec  =  mequtorv(p_o,  f_o,  g_o,  h_o,  k_o,  L_o) 

7.  This  function  inputs  the  initial  modified  equinoctial  element  set  (p,  f, 

7.  g,  h,  k,  &  L)  and  outputs  the  initial  state  vector  x_o  (composed  of  the  r 
7.  and  v  vectors)  based  on  the  inital  element  set. 

global  g_s  G  M_earth  R_earth  mu_earth 

alpha  =  sqrt(h_o~2  -  k_o~2) ;  s  =  sqrt(l  +  h_o~2  +  k_o~2) ;  w  =  1  + 
f_o*cos(L_o)  +  g_o*sin(L_o) ;  r  =  p_o/w; 

xo  =  (r/s~2) * (cos (L_o)  +  alpha~2*cos (L_o)  +  2*h_o*k_o*sin(L_o) ) ; 
yo  =  (r/s~2) * (sin(L_o)  -  alpha~2*sin(L_o)  +  2*h_o*k_o*cos (L_o) ) ; 
zo  =  (2*r/s~2) * (h_o*sin(L_o)  -  k_o*cos (L_o) ) ; 

xdoto  =  (-l/s~2) *sqrt (mu_earth/p_o) * (sin(L_o)  + 
alpha~2*sin(L_o) . . . 

-  2*h_o*k_o*cos (L_o)  +  g_o  -  2*f _o*h_o*k_o  +  alpha~2*g_o) ; 
ydoto  =  (-l/s~2) *sqrt (mu_earth/p_o) * (-cos (L_o)  + 
alpha~2*cos (L_o) . . . 

+  2*h_o*k_o*sin(L_o)  -  f_o  +  2*g_o*h_o*k_o  +  alpha~2*f_o) ; 
zdoto  =  (2/s~2) *sqrt (mu_earth/p_o) * (h_o*cos (L_o)  +  k_o*sin(L_o) . .  . 

+  f _o*h_o  +  g_o*k_o) ; 

x_o_vec  =  [xo;  yo;  zo;  xdoto;  ydoto;  zdoto]; 
return 


C.9  oblateearth.m 

function  J2term  =  oblateearth(t ,x) 

global  SF  c  g_s  G  M_earth  R_earth  mu_earth  m  epoch_yr  A_plate 
A_cyl . . . 

Earth_aphelion_table  c_R_plate  c_R_cyl  J2 
e_earth  =  sqrt (0 . 006694385) ; 
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p  =  x(l) ;  f  =  x(2) ;  g  =  x(3);  h  =  x(4);  k  =  x(5);  L  =  x(6); 


7.  compute  r  vector  from  EQU  set 

alpha  =  sqrt(h~2  -  k~2)  ;  s  =  sqrt(l  +  h~2  +  k~2)  ;  w  =  1  +  f*cos(L) 
+  g*sin(L) ;  r  =  p/w;  ri  =  (r/s~2) * (cos (L)  +  alpha~2*cos (L)  + 
2*h*k*sin(L) ) ;  rj  =  (r/s~2) * (sin(L)  -  alpha~2*sin(L)  + 

2*h*k*cos (L) ) ;  rk  =  (2*r/s~2)*(h*sin(L)  -  k*cos(L));  r.vec  =  [ri; 
r  j  ;  rk]  ; 

xdot  =  (-l/s~2)*sqrt(mu_earth/p)*(sin(L)  +  alpha~2*sin(L) . . . 

-  2*h*k*cos(L)  +  g  -  2*f*h*k  +  alpha~2*g) ; 
ydot  =  (-l/s~2) *sqrt (mu_earth/p) * (-cos (L)  +  alpha~2*cos (L) . . . 

+  2*h*k*sin(L)  -  f  +  2*g*h*k  +  alpha~2*f); 
zdot  =  (2/s~2) *sqrt (mu_earth/p) * (h*cos (L)  +  k*sin(L)  +  f*h  +  g*k) ; 
v_vec  =  [xdot;  ydot;  zdot]; 

ir  =  r_vec/norm(r_vec) ;  in  = 

cross (r_vec , v_vec) /norm(cross (r_vec , v_vec) ) ;  it  =  cross(in, ir)  ; 

R  =  [ir  it  in] ;  %  Rotation  matrix  from  ECI  frame  to  RTN  frame 

a  =  R_earth; 


%  find  the  geocentric  latitude  (Vallado  pg  146) 

phi_gc  =  atan2  (rk,  sqrt  (ri~2  +  rj~2));  7.  rad,  pg  145  instead 

7,  define  the  Legendre  polynomials  (Vallado  pg  517) 

P2  =  (l/2)*(3*(sin(phi_gc))~2  -  1);  7.  P_2,0 
P2prime  =  3*sin(phi_gc) ; 

7.  Dr.  Wiesel’s  solution,  interpreted  from  Betts 


7.  find  the  perturbing  acceleration 

delta_gr  =  (-mu_earth/r~2)*(3*(a/r) ~2*P2*J2) ;  delta_gn  = 
(-mu_earth*cos (phi_gc) /r~2) * ( (a/r) ~2*P2prime*J2) ; 

en_vec  =  [001]’;  inorth  = 

(en_vec- (en_vec ’ *ir) *ir) /norm(en_vec-(en_vec ’ *ir) *ir) ; 

7.  describes  the  local  North  direction 

delta_g  =  delta_gn*inorth  -  delta_gr*ir; 

7.  gravitational  disturbing  acceleration 

J2term  =  R’*delta_g; 
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return 


C.10  propxequ.m 

function  xdot  =  propxequ(t ,x) 

V.  This  function  inputs  the  initial  state  vector  x_o_vec  (composed  of  r_o 
%  and  v_o  vectors)  and  outputs  the  state  vector  x_vec  (composed  of  r  and 
1  vectors)  based  on  those  parameters  and  applicable  perturbations. 

global  n  g_s  G  M_earth  R_earth  mu_earth  rho  J2  SF  c  m  epoch_yr 
epoch_day. . . 

epoch_frac  aa  countt 

%  track  progress : 

countt=countt+l ;  orbitno  =  t/9.037752480834280e+004;  if  orbitno  >= 
aa 

aa;  °/,  uncomment  to  track  compile  progress 
aa  =  aa+1; 

end 

p  =  x(l) ;  f  =  x(2) ;  g  =  x(3);  h  =  x(4);  k  =  x(5);  L  =  x(6); 


°/0  calculate  current  r_vec,  v_vec  and  R 

alpha_squared  =  h~2  -  k~2;  w  =  1  +  f*cos(L)  +  g*sin(L) ;  r  =  p/w; 

Chi  =  sqrt(h"2  +  k~2) ;  s_squared  =  1  +  Chi~2; 

ri  =  (r/s_squared) * (cos (L)  +  alpha_squared*cos (L)  +  2*h*k*sin(L) ) ; 
rj  =  (r/s_squared) * (sin(L)  -  alpha_squared*sin(L)  +  2*h*k*cos (L) ) ; 
rk  =  (2*r/s_squared)*(h*sin(L)  -  k*cos(L)); 

r_vec  =  [ri;  r j ;  rk] ;  %  ECI  position  vector  of  the  satellite  (km) 

xdot  =  (-l/s_squared) *sqrt (mu_earth/p) * (sin(L)  + 
alpha_squared*sin(L) . . . 

-  2*h*k*cos (L)  +  g  -  2*f*h*k  +  alpha_squared*g) ; 
ydot  =  (-l/s_squared) *sqrt (mu_earth/p) * (-cos (L)  + 
alpha_squared*cos (L) . . . 

+  2*h*k*sin(L)  -  f  +  2*g*h*k  +  alpha_squared*f ) ; 
zdot  =  (2/s_squared) *sqrt (mu_earth/p) * (h*cos (L)  +  k*sin(L)  +  f*h  + 
g*k) ; 

v_vec  =  [xdot;  ydot;  zdot];  °/0  ECI  velocity  vector 
x_vec  =  [r_vec;  v_vec] ; 

ir  =  r_vec/norm(r_vec) ;  in  = 

cross (r_vec , v_vec) /norm(cross (r_vec , v_vec) ) ;  it  =  cross(in, ir) ; 

R  =  [ir  it  in]  ;  °/0  Rotation  matrix  from  ECI  frame  to  RTN  frame 


%  calculate  current  Julian  date  &  Earth-Sun  vector 
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JD  =  juldatet(t); 


7.  JD  is  in  ref  to  UTC  (let  JD_UTC  =  JD_UT1) 


r_ES_vec  =  earthsunt (t) ;  %  ECI  position  vector  of  the  Sun  (km) 


%  calculate  matrix  form  A  &  b  (3  Nov  05  changes) 

A  =  [0  (2*p/w) *sqrt (p/mu_earth)  0;... 

sqrt (p/mu_earth) *sin(L)  sqrt(p/mu_earth)*(l/w)*((w+l)*cos(L)+f)  . . 
sqrt(p/mu_earth)*(g/w)*(h*sin(L)-k*cos(L)) ; . . . 

-sqrt (p/mu_earth) *cos (L)  sqrt(p/mu_earth)*(l/w)*((w+l)*sin(L)+g) . . . 
sqrt(p/mu_earth)*(f/w)*(h*sin(L)-k*cos(L)) ;  .  .  . 

0  0  sqrt (p/mu_earth) * (s_squared*cos (L) / (2*w) ) ; . . . 

0  0  sqrt (p/mu_earth) * (s_squared*sin(L) / (2*w) ) ; . . . 

0  0  sqrt(p/mu_earth)*(l/w)*(h*sin(L)-k*cos(L))] ; 

b=  [0;  0;  0;  0;  0;  sqrt (mu_earth*p) * (w/p) ~2] ; 


1  J2  effects 

J2term  =  oblateearth(t ,x) ;  %  same  results  as  Vallado! 

J2_r  =  J2term(l,:);  °/«  radial  component  of  SRP  perturbing  acceleration  vector 

J2_t  =  J2term(2,:);  °L  tangential  component  of  SRP  perturbing  acceleration  vector 

J2_n  =  J2term(3,:);  °/«  normal  component  of  SRP  perturbing  acceleration  vector 


%  Solar  third-body  effects 

mu_sun  =  1 . 32712428ell ;  °/«  km~3/s~2 

da  =  - (mu_sun/norm(r_ES_vec) ~3) * (r_vec  -  ... 

3*r_ES_vec*(dot(r_vec,r_ES_vec)/norm(r_ES_vec)~2)  -  ... 
15*r_ES_vec* (dot (r_vec ,r_ES_vec) /norm(r_ES_vec) ~2) ~2) ; 
a_sun  =  R’*da;  sun_r  =  a_sun(l,:);  sun_t  =  a_sun(2,:);  sun_n  = 
a_sun(3, : ) ; 


%  find  the  number  of  Julian  centuries 

T_TDB  =  (JD  -  2451545.0)736525;  7,  Vallado  pg  188 


7,  SRP  effects  (only  in  sat-sun  radial  direction  according  to  Vallado) 

7.  implement  shadow  algorithm  (page  297  of  Vallado) 
tau_min  =  (norm(r_vec) ~2-dot (r_vec ,r_ES_vec) ) / . . . 

(norm(r.vec) ~2+norm(r_ES_vec) ~2-2*dot (r_vec ,r_ES_vec) ) ; 

7«  make  sure  vectors  have  ssime  units! 

min_val2  =  ( (l-tau_min) *norm(r_vec) ~2  +  r_vec’*r_ES_vec*tau_min)/6378~2; 
7.  must  be  in  ER~2! 

if  ((tau_min  <  0.0)  I  (tau_min  >1)) 
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%  satellite  is  illuminated;  run  SRP  algorithm 
SRP  =  srpaccel(JD,x_vec,R,r_ES_vec) ; 

SRP_r  =  SRP ( 1 , : ) ; 

7.  radial  component  of  SRP  perturbing  acceleration  vector 
SRP_t  =  SRP (2, :) ; 

7«  tangential  component  of  SRP  perturbing  acceleration  vector 
SRP_n  =  SRP (3, :)  ; 

7.  normal  component  of  SRP  perturbing  acceleration  vector 
shadow  =  1 ; 

elseif  min_val2  >=  1.0 

7.  satellite  is  illuminated;  run  SRP  algorithm 
SRP  =  srpaccel(JD,x_vec,R,r_ES_vec) ; 

SRP_r  =  SRP ( 1 , : ) ; 

7,  radial  component  of  SRP  perturbing  acceleration  vector 
SRP_t  =  SRP (2, :) ; 

7«  tangential  component  of  SRP  perturbing  acceleration  vector 
SRP_n  =  SRP (3, :) ; 

7.  normal  component  of  SRP  perturbing  acceleration  vector 
shadow  =  1 ; 

else 

7.  satellite  is  NOT  illuminated;  skip  SRP  algorithm  since  not  affected 
SRP_r  =  0; 

7o  radial  component  of  SRP  perturbing  acceleration  vector 
SRP_t  =  0; 

7«  tangential  component  of  SRP  perturbing  acceleration  vector 
SRP_n  =  0 ; 

7.  normal  component  of  SRP  perturbing  acceleration  vector 
shadow  =  0 ; 

end 


7,  formulate  perturbation  acceleration  vector  (RTN  frame) 

Pr  =  J2_r  +  SRP_r  +  sun_r;  Pt  =  J2_t  +  SRP_t  +  sun_t ;  Pn  =  J2_n  + 
SRP_n  +  sun_n; 

P  =  [Pr;  Pt;  Pn] ; 

7,  state  equation  solution 

xdot  =  A*P  +  b; 

return; 


C.ll  srpaccel.  m 

function  SRP  =  srpaccel(JD,x_vec,R,r_ES_vec) 

7.  This  function  inputs  the  iteration  time  step  and  state  vector  and  output 
7.  the  total  disturbing  acceleration  caused  by  SRP  forces, 
global  SF  c  g_s  G  M_earth  R_earth  mu_earth  m  epoch_yr  A_plate 
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A_cyl . . . 

Earth_aphelion_table  c_R_plate  c_R_cyl  cont 

p_SR_simple  =  SF/c;  %  simple  approximation  (W/km~2)/ (km/sec) =W*km/sec 
row  =  epoch_yr-2004;  D_aph  =  Earth_aphelion_table (row, 2) ; 
p_SR_complex  =  (1358/ (1 . 004+0 . 0334*cos (D_aph) )* (1/c) ) ; 

7.  more  complex  approximation  (page  547) 

C_Rb  =  2;  7.  coefficient  of  reflectivity  of  s/c  body  (dimensionless) 

7.  2  indicates  flat  mirror  perp  to  sun 
%  1  indicates  a  black  body  surface 
C_Rp  =  2; 

7.  coefficient  of  reflectivity  of  parabolic  reflectors  (dimensionless) 

r_vec  =  [x_vec(l);  x_vec(2);  x_vec(3)]; 

7.  this  is  the  Earth-Satellite  vector 
v_vec  =  [x_vec(4);  x_vec(5);  x_vec(6)]; 

r_satsun  =  -r_vec  +  r_ES_vec;  7,  -r_earthsat  +  r_earthsun 

switch  cont 

case  1  °l  controller  1 


7,  Baseline  case  (no  controller,  area  of  despun  platform  =  zero) 
p_SR_simple  =  SF/c; 

7.  simple  approximation  (W/km~2) /(km/sec)=W*km/sec 
row  =  epoch_yr-2004; 

D_aph  =  Earth_aphelion_table(row,2) ; 

p_SR_complex  =  (1358/ (1 . 004+0 . 0334*cos (D_aph) ) *(l/c) ) ; 

7.  more  complex  approximation  (page  547) 

C_Rb  =  1.5; 

7.  coefficient  of  reflectivity  of  s/c  body  (dimensionless) 

7.  2  indicates  flat  mirror  perp  to  sun 
7.  1  indicates  a  black  body  surface 

ir  =  r_vec/norm(r_vec) ; 

in  =  cross (r_vec , v_vec) /norm(cross (r.vec , v_vec) ) ; 
it  =  cross (in, ir) ; 

R  =  [ir  ’  ;  it  ’  ;  in’ ] ; 

7.  Rotation  matrix  from  ECI  frame  to  RTN  frame 

r_satsun  =  -r_vec  +  r_ES_vec;  7,  -r_earthsat  +  r_earthsun 
s_vec  =  r_satsun/norm(r_satsun)  ;  °L  sat-sun  UNIT  vector 

7.  for  the  rotating  cylinder: 
phi_cyl  =  0; 

7.  reflected  or  incident  angle  of  the  cylinder 
7.  (always  perpendicular  to  s_vec  since  a  cylinder) 
n_vec_c  =  s_vec; 

7.  again,  face  of  cylinder  always  perpendicular  to  s_vec 
n_vec_cyl  =  n_vec_c/norm(n_vec_c) ;  %  UNIT  vector 
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I  for  the  despun  platform: 

A_expp  =  0; 
phi_plate  =  0; 
n_vec_plate  =  n_vec_cyl; 

J  note:  reflectivity  coefficients  [absorption;  diffuse;  specular] 
a_SRP_plate  =  (p_SR_complex*A_expp*cos(phi_plate)/m)* . . . 

(2* ( (c_R_plate (2) /3)  +  c_R_plate (3) *cos (phi_plate) ) * . . . 
n_vec_plate  +  (l-c_R_plate(3))*s_vec) ; 
a_SRP_cyl  =  (p_SR_complex*A_cyl*cos(phi_cyl)/m)* . . . 

(2*((c_R_cyl(2)/3)  +  c_R_cyl (3) *cos (phi_cyl) ) *n_vec_cyl  +... 
(l-c_R_cyl(3))*s_vec) ; 
a_SRP  =  -a_SRP_plate  -  a_SRP_cyl; 

SRP  =  R*a_SRP;  7«  SRP  acceleration  vector  in  RTN  frame 

case  2  7.  controller  2 

*/,  Vallado  (page  549) 

s_vec  =  r_satsun/norm(r_satsun) ;  7«  sat-sun  UNIT  vector 

norm(s_vec) ; 

costheta  =  dot (r_satsun, v_vec) . /(norm(r_satsun)*norm(v_vec)) ; 

7.  simple  on/off  control  algorithm 

if  costheta  <=  0  7o  full  plate  area  toward  sun 

gamma  =  acos  (dot  (r_vec ,  s_vec) /norm(r_vec) )  ;  7o  radians 
A_expp  =  A_plate;  °L  km~2,  flat  plate  approximation  3m  x  1.1176m 
else  7.  no  plate  area  toward  sun 

gamma  =  acos  (dot  (r_vec ,  s_vec) /norm(r_vec)  )  -  pi/2;  7o  radians 
A_expp  =  0;  7o  km~2,  flat  plate  approximation 

end 

theta_vallado  =  acos(dot(r_vec,v_vec)/ (norm(r_vec) *norm(v_vec) ) ) ; 

7.  radians 

7.  for  the  antennae: 

phi_plate  =  acos (dot (r_vec , s_vec) /norm(r_vec) ) -gamma; 

7.  reflected  or  incident  angle  of  the  antennae 


7.  calculate  normal  vector,  n_vec  (based  on  commanded  angle) 

7.  calculate  n_vec  (rotate  from  r_vec  to  n_vec  through  gamma)  : 
R_gamma  =  [cos (gamma)  sin (gamma)  0; -sin (gamma)  cos (gamma)  0;  0  0  1] 
n_vec_p  =  R’ *R_gamma*R*r_vec ; 

n_vec_plate  =  n_vec_p/norm(n_vec_p) ; 

7.  UNIT  vector  perpendicular  to  plate  surface,  directed  opposite 
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7.  of  normal  force 


(thus  towards  Sun) 


case  3  7.  controller  3 


7.  Coverstone  model  for  control  algorithm  (modified) 

s_vec  =  r_satsun/norm(r_satsun) ; 

7.  sat-sun  UNIT  vector  (km)  in  ECI  frame 

7.  1.  calculate  the  number  of  Julian  centuries 
T_UT1  =  (JD  -  2451545. 0)/36525;  */,  Vallado  pg  188 
lambda_M_sundeg  =  280.460  +  36000 . 770*T_UT1 ; 

*/,  mean  ecliptic  longitude  of  the  Sun  (deg)  ,  Vallado  pg  265 
T_TDB  =  T_UT1 ; 

M_sundeg  =  35T.52TT233  +  35999 . 05034*T_TDB; 

*/,  Sun’s  mean  anomaly  (deg) 

lambda_ecliptic  =  lambda_M_sundeg  +  1 . 914666471*sind(M_sundeg)  +.. 
0 . 019994643*sind(2*M_sundeg) ; 

7.  ecliptic  longitude  of  the  Sun  (deg) 

•/,  2.  rotate  v_vec  from  ECI  to  heliocentric  XYZ  frame: 

R_XYZf romECI  =[10  0;0  cosd(23.5)  sind(23 . 5) ; .  .  . 

0  -sind(23.5)  cosd(23.5)]; 
v_vec_XYZ  =  R_XYZfromECI*v_vec; 

•/,  3.  rotate  v_vec_XYZ  to  alpha-beta-z  (ABC)  frame: 

R_ABZfromXYZ  =  [cosd(lambda_ecliptic)  sind(lambda_ecliptic)  0; 

-sind(lambda_ecliptic)  cosd(lambda_ecliptic)  0;  001]; 
v_vec_ABZ  =  R_ABZfromXYZ*v_vec_XYZ ; 

7.  4.  calculate  the  desired  normal  vector,  n_vec_ABZ: 
zeta  =  (-3*v_vec_ABZ(l)*v_vec_ABZ(2)  -  v_vec_ABZ(2) * . . . 
sqrt (9*v_vec_ABZ(l) ~2  +  8*(v_vec_ABZ(2)~2  +... 
v_vec_ABZ(3) ~2) ) ) / (4* (v_vec_ABZ(2) ~2  +  v_vec_ABZ(3) ~2) ) ; 
n_alpha  =  -abs (v_vec_ABZ(2) ) /sqrt (v_vec_ABZ(2) ~2  +... 

zeta~2*(v_vec_ABZ(2)~2  +  v_vec_ABZ(3) ~2)) ; 
n_beta  =  zeta*n_alpha; 

n_z  =  (v_vec_ABZ(3)/v_vec_ABZ(2))*n_beta; 
n_vec_ABZ  =  -[n_alpha;  n_beta;  n_z] ; 

7.  5.  rotate  n_vec_ABZ  from  ABZ  frame  to  heliocentric  XYZ  frame: 
n_vec_XYZ  =  R_ABZfromXYZ’ *n_vec_ABZ; 

7.  6.  rotate  n_vec_XYZ  from  XYZ  frame  to  ECI  frame: 

n_vec  =  R_XYZf romECI ’ *n_vec_XYZ;  %  unit  normal  vector  in  ECI  frame 

n_vec_plate  =  n_vec; 

phi_plate  =  acos (dot (n_vec , s_vec) ) ; 

7.  reflected  or  incident  angle  of  the  antennae  (rad) 
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A_expp  =  A_plate; 

*/.  assume  controller  faces  full  area  of  plate  toward  Sun  always 

end 

switch  cont 

case  {2,3}  70controller  2  or  3 
7.  for  the  rotating  cylinder: 
phi_cyl  =  0; 

7.  reflected  or  incident  angle  of  the  cylinder 
7«  (always  perpendicular  to  s_vec  since  a  cylinder) 
n_vec_c  =  s_vec; 

7.  again,  face  of  cylinder  always  perpendicular  to  s_vec 
n_vec_cyl  =  n_vec_c/norm(n_vec_c) ;  %  UNIT  vector 

7.  note:  reflectivity  coefficients  [absorption;  diffuse;  specular] 
a_SRP_plate  =  (p_SR_complex*A_expp*cos(phi_plate)/m)* . . . 

(2* ( (c_R_plate (2) /3)  +  c_R_plate (3) *cos (phi_plate) ) * . . . 
n_vec_plate  +  (l-c_R_plate(3))*s_vec) ; 
a_SRP_cyl  =  (p_SR_complex*A_cyl*cos (phi_cyl) /m) * (2* ( (c_R_cyl (2) /3) 
+  c_R_cyl (3) *cos (phi_cyl) ) *n_vec_cyl  +  (l-c_R_cyl(3))*s_vec) ; 
a_SRP  =  -a_SRP_plate  -  a_SRP_cyl; 

SRP  =  R’*a_SRP;  %  SRP  acceleration  vector  in  RTN  frame 

end 

return; 
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Appendix  D.  Disposal  Orbit  Guidelines 


D.l  Satellite  Disposal  Guidelines 

The  United  States  Space  Command  (USSPACECOM)  is  tasked  with  defining  and  directing 
policies  with  respect  to  the  disposal  of  satellites  owned  by  the  DoD.  The  pertaining  policy  directive, 
UPD10-39  (reference  Section  D.2),  states  that  it  is  critical  to  remove  abandoned  spacecraft  from 
orbit  at  the  end  of  their  orbital  life.  In  “high  value”  operational  orbits,  there  is  an  increasing 
threat  of  collision  between  satellites.  Removing  unused  satellites  from  a  given  orbit  (particularly 
low-Earth  orbit  and  geosynchronous  orbit)  would  both  reduce  the  number  of  potential  collisions 
and  allow  additional  room  for  new  satellites  to  be  launched. 

Necessary  paperwork  and  approvals  aside,  the  portion  of  the  document  pertinent  to  this  thesis 
lies  in  Paragraph  5.4,  entitled  “Above  Geosynchronous  Orbit”.  Although  the  document  states  that 
a  satellite  may  be  disposed  of  using  any  of  the  methods  listed  in  Section  5,  the  current  orbit  of  the 
DSCS  II  satellite  leads  toward  implementation  of  Paragraph  5.4.  To  summarize,  USSPACECOM 
directs  that  the  abandoned  satellite  must  be  maneuvered  to  an  orbit  with  a  perigee  altitude  above 
36100  km.  The  DSCS  II  F-13  satellite  already  exceeds  this  requirement,  however,  for  the  purposes 
of  this  thesis  it  will  still  be  used  to  determine  how  long  it  would  take  to  raise  its  orbit  a  specified 
amount. 

D.2  UPD 10-39 
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BY  ORDER  OF  USCINC SPACE 


UNITED  STATES  SPACE  COMMAND 
POLICY  DIRECTIVE  10-39 

1  FEBRUARY  2001 
Operations 

SATELLITE  DISPOSAL  PROCEDURES 


NOTICE:  This  publication  is  available  digitally  at:  https://midway.peterson.af.mil/21etters/sc/css/scr/ 
norad-us/utable .  htm . 


OPR:  J35X  (LCDR  Steven  R.  White)  Certified  by:  JS  (COL  Robert  A.  Hammerle) 

Supersedes  UPD10-39,  3  Nov  97.  Pages:  7 

Distribution:  F 


This  policy  directive  establishes  United  States  Space  Command  (USSPACECOM)  policy  on  the  disposal 
of  Department  of  Defense  owned  satellites  which  Commander-in-Chief,  US  Space  Command  (USCINC- 
SPACE)  exercises  Combatant  Command  (COCOM)  authority  as  defined  by  the  current  Forces  for  Uni¬ 
fied  Commands  document.  It  does  not  apply  to  Air  Force  Reserve  Command  (AFRC)  nor  Air  National 
Guard  (ANG)  units. 

SUMMARY  OF  REVISIONS 

This  policy  directive  incorporates  updated  guidance  for  satellite  disposal  methods/regions,  with  majority 
of  revisions  found  in  paragraphs  5.-5. 6.  A  bar  (I)  indicates  a  revision  from  the  previous  edition. 


1.  General.  The  space  environment  is  critical  as  a  medium  for  the  collection  and  dissemination  of  data 
and  information  crucial  to  the  warfighter’ s  mission  accomplishment.  If  the  practice  of  abandoning  space¬ 
craft  in-place  within  high  value  operational  orbits  at  the  end  of  their  mission  life  continued,  there  will  be 
an  ever-increasing  threat  of  collision  between  abandoned  and  operational  satellites.  Furthermore,  such 
collisions  themselves  will  become  a  major  secondary  source  of  orbital  debris,  posing  a  geometrically 
growing  threat  to  future  space  operations  that  will  be  virtually  impossible  to  control.  To  protect  this  vital 
resource  from  being  polluted  by  space  debris,  this  document  provides  direction  and  guidance  for  the 
proper  disposal  of  satellites.  These  procedures  comply  with  the  minimization  and  mitigation  of  space 
debris  as  directed  by  National  and  Department  of  Defense  (DOD)  Space  Policy.  The  guidelines  do  not 
preclude  any  end-of-life  testing  that  the  organization  with  satellite  system  responsibility  deems  necessary 
either  prior  to  or  following  the  placement  into  a  disposal  orbit. 

2.  Mission  Payload/Vehicle  Health.  Removing  a  non-mission  capable  satellite  from  its  operational 
orbit  into  an  established  disposal  region  is  of  paramount  importance.  As  a  satellite  approaches  the  end  of 
its  operational  life,  each  SATCOM  System  Expert  (SSE)  responsible  for  the  satellite  bus,  or  their  equiva¬ 
lent  for  non-SATCOM  systems  (see  Table  1.  and  Table  2.),  will  ensure  that  every  satellite  maintains  ade¬ 
quate  disposal  capability.  This  is  to  include  assuring  command/control  capability  and  maintaining  the 
required  amount  of  fuel  to  reach  the  disposal  region.  Disposal  of  vehicles  approaching  the  end  of  their 
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operational  life  will  be  recommended  if  further  degradation  precludes  future  removal  from  high-value 
operational  orbits. 


Table  1.  Organizations  With  System  Responsibility  (Communication  Satellites  and  Secondary 
Payloads). 


SATELLITE 

SATCOM 

SATCOM 

SYSTEMS/ 

OPERATIONAL 

SYSTEM 

PAYLOADS 

MANAGER  (SOM) 

EXPERT  (SSE) 

FLTSAT,  UHF  Follow-On  (UFO), 
FEP,  UFO(E),  UFO(EE),  Polar  EHF, 
GBS  Phase  11 

USSPACECOM 

NA  V  SPACECOM 

DSCS 

USSPACECOM 

DISA 

Milstar,  AFSATCOM.  SCTS,  Polar 
UHF,  LES 

USSPACECOM 

AFSPC 

Table  2.  Organizations  With  System  Responsibility  (Non-Communication  Satellites). 


SATELLITE 

SYSTEMS 

ORGANIZATION 

DSP 

SPACEAF 

DMSP 

SPACEAF  in  coordination  with  the  National  Polar-Orbiting  Opera¬ 
tional  Environment  Satellite  System  (NPOESS)  Integrated  Program 
Office  (IPO) 

GPS 

SPACEAF 

3.  Disposal  Criteria.  SATCOM  System  Experts  (SSE's),  or  their  equivalent  for  non-SATCOM  systems, 
will  coordinate  with  sustainment  and  satellite  operations  activities,  components  or  agencies  with  second¬ 
ary  payloads  on  the  bus,  and  develop  specific  criteria  for  satellite  disposal.  USSPACECOM  components 
and  agencies  will  forward  the  criteria  to  USSPACECOM/J3  (and  J6  for  SATCOM  systems)  for  review 
and  approval.  The  disposal  criteria  will  include  the  minimum  acceptable  levels  of  bus  support  to  the  pay- 
loads,  payload  capability  and  capacity  (primary  and  secondary),  vehicle  command/control  capability, 
vehicle  power  capacity,  fuel  requirements  for  disposal  maneuver,  and  any  other  disposal  maneuver 
requirements. 

4.  Standard  Operations.  As  part  of  standard  operations,  SATCOM  System  Experts  (SSE's),  or  their 
equivalent,  will  monitor  satellite  capability  criteria.  Once  a  satellite  has  been  designated  non-mission 
capable  per  the  established  criteria  and/or  is  demonstrating  potential  disposal  capability  problems,  the 
organization  will  forward  a  disposal  recommendation  as  described  below.  At  a  minimum,  the  recommen¬ 
dation  will  include  the  disposal  criteria  that  are  being  met  and  projected  date  of  disposal.  Standard  dis¬ 
posal  recommendations  should  be  sent  three  months  in  advance  of  the  anticipated  disposal  date. 

4.1.  Military  Satellite  Communications  (MILSATCOM)  Disposal  Recommen  dation.  SAT¬ 
COM  System  Experts  (SSE’s)  will  forward  the  satellite  disposal  recommendation  to  USSPACECOM/ 

J6.  USSPACECOM/J3/J6  will  then  coordinate  with  other  government  agencies  to  determine  any 
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requirements  for,  and  the  feasibility  of,  potential  alternate  use.  USSPACECOM  will  then  consolidate 
these  requirements,  and  make  a  recommendation  to  Joint  Staff/J6  on  the  satellite’s  disposition.  Once 
final  approval  is  received,  USSPACECOM/J3  will  direct  disposal  action.  The  SSE  will  develop  and 
provide  a  plan  of  action  and  timeline  for  disposal/transfer  to  include  end-of-life  testing  requirements, 
Test  and  Checkout  (TACO),  and  disposal  orbit  parameters  to  USSPACECOM/J6/Global  SATCOM 
Support  Center  (GSSC).  A  final  status  report  will  be  released  by  the  SATCOM  Operational  Manager 
(SOM)  to  the  SATCOM  community  (i.e.,  CINCs,  Joint  Program  Offices,  USSPACECOM  compo¬ 
nents,  etc.).  If  the  satellite  is  recommended  for  continued  use  by  another  agency,  then  the  operations/ 
maintenance/support  and  final  disposal  plan  will  be  coordinated  with  the  Joint  Staff  J3/J6. 

4.2.  Non-MILSATCOM  Disposal  Recommendation.  The  organization  with  satellite  system 
responsibility  will  submit  their  disposal  recommendation  to  USSPACECOM/J3,  who  will  make  the 
determination  for  final  disposition.  Once  approval  has  been  granted,  the  responsible  organization  will 
develop  and  provide  a  plan  of  action  and  timeline  for  disposal  to  include  end-of-life  testing  require¬ 
ments,  TACO,  and  disposal  orbit  parameters  to  USSPACECOM/J3. 

4.3.  Emergency  Situation.  In  the  case  of  an  emergency  where  it  is  believed  that  there  is  less  than  30 
days  to  react,  the  decision  process  for  the  safe  disposal  of  a  satellite  must  be  expedited.  The  Satellite 
Control  Authority  (SCA)/system  operator  will  make  a  telephone  request,  followed  by  message,  for 
disposal  to  the  SATCOM  System  Experts  (SSE),  or  their  equivalent  for  non-SATCOM  systems.  The 
SATCOM  System  Experts  (SSE),  or  their  equivalent,  will  coordinate  via  phone  and  message  with 
USSPACECOM/J6  for  MILSATCOM  systems  or  USSPACECOM/J3  for  non-MILSATCOM  sys¬ 
tems.  USSPACECOM/J3/J6  will  coordinate  the  disposal  recommendation. 

4.4.  Maneuver  Vector  Screening.  The  component  will  ensure  planned  disposal  vectors  and  orbit 
parameters  are  submitted  to  Cheyenne  Mountain  Operations  Center  (CMOC),  J3SY  Combat  Analysis 
Branch  before  disposal  for  approval  of  reentry  locations  and/or  orbital  safety  screening  for  possible 
conflicts.  USSPACECOM/J3  will  resolve  any  conflicts  concerning  reentry  locations  and/or  orbital 
safety  screening.  After  disposal  maneuvers  have  been  completed,  the  component  will  ensure  that  the 
final  vectors  are  provided  to  CMOC.  If  the  satellite  fails  to  reach  proper  disposal  orbit,  the  Space 
Control  Center  (SCC)  will  conduct  conjunction  analysis.  If  there  is  a  possible  conflict  with  a  manned 
space  vehicle  or  other  objects,  continue  routine  conjunction  analysis.  Notify  the  Space  Operations 
Center  (SPOC)  and  other  space  partners  of  all  conjunctions  associated  with  disposal  operations. 

5.  Post  Mission  Disposal.  Satellites  should  be  disposed  of  by  one  of  the  following  methods.  Because  of 
fuel  gauging  uncertainties  near  the  end  of  mission,  a  disposal  plan  should  use  a  maneuver  strategy  that 
most  reduces  the  risk  of  leaving  the  satellite  near  a  high-value  orbit.  The  disposal  region  perigee  and  apo¬ 
gee  altitudes  listed  below  should  prevent  lunar  and  solar  perturbations  from  causing  these  satellites  to 
interfere  with  high-value  operational  orbits.  Programs  whose  current  design  capabilities  are  unable  to 
meet  the  guidelines  listed  in  paragraphs  5.1. -5.6.,  will  optimize  disposal  actions  based  on  the  published 
design  capabilities. 

5.1.  Atmospheric  Reentry.  Maneuver  the  satellite  to  an  orbit  in  which,  using  conservative  projec¬ 
tions  for  solar  activity  and  other  perturbations,  atmospheric  drag  will  cause  atmospheric  reentry 
within  25  years  after  completion  of  mission.  If  atmospheric  reentry  is  performed  by  a  planned  deorbit, 
it  will  be  planned  such  that  any  remaining  portions  of  the  satellite  will  impact  the  earth  only  in 
non-populated,  preferably  oceanic,  areas. 
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5.2.  Between  Low  Earth  Orbit  (LEO)  and  Medium  Earth  Orbit  (MEO)/Semisynchronous. 

Maneuver  the  satellite  to  an  orbit  with  perigee  altitude  above  2,000  km  and  apogee  altitude  below 
19,700  km. 

5.3.  Between  Medium  Earth  Orbit  (MEO)/Semisynchronous  and  Geosynchronous  Earth  Orbit 
(GEO).  Maneuver  the  satellite  to  an  orbit  with  perigee  altitude  above  20,700  km  and  apogee  altitude 
below  35,300  km. 

5.4.  Above  Geosynchronous  Earth  Orbit  (GEO).  Maneuver  the  satellite  to  an  orbit  with  perigee 
altitude  above  36,100  km. 

5.5.  Heliocentric,  Earth-escape.  Maneuver  to  remove  the  satellite  from  Earth  orbit,  into  a  heliocen¬ 
tric  orbit. 

5.6.  Direct  retrieval.  Retrieve  the  satellite  and  remove  it  from  orbit  as  soon  as  practical  after  mission 
completion. 

6.  Safing.  Properly  safing  the  payloads  and  the  bus  is  a  critical  step  in  the  disposal  process.  This  process 
will  assist  in  the  elimination  of  stored  energy  from  the  satellite  and  limit  the  probability  of  post  mission 
explosion  spreading  more  debris.  Safing  procedures  to  be  considered  vary  with  each  system,  but  may 
include  burning  all  residual  fuel  to  depletion  and  leaving  fuel  lines  with  valves  open;  disabling  all  battery 
charging  systems;  leaving  batteries  in  a  permanent  discharge  state;  venting  all  pressurized  systems; 
removing  power  from  control  moment  gyroscopes;  disabling  transmitters;  deactivating  range  safety  sys¬ 
tems;  depleting  all  volatile  liquids;  or  stabilizing  the  spacecraft  in  a  neutral  thermal  flight  mode. 

7.  Review.  This  policy  directive  will  be  reviewed  at  least  every  3  years  to  ensure  the  most  current  guid¬ 
ance  is  provided. 


RALPH  E.  EBERHART,  General,  USAF 
Commander  in  Chief 
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Attachment  1 

GLOSSARY  OF  REFERENCES  AND  SUPPORTING  INFORMATION 


References 

Presidential  Decision  Directive  (PDD)  NSC-49/NSTC-8,  National  Space  Policy 
Department  of  Defense  Directive  3100.10,  Space  Policy 

United  States  Space  Command  Instruction  (UI)  13-4,  Minimization  and  Mitigation  of  Space  Debris 
NASA  Safety  Standard  1740 

National  Reconnaissance  Office  Instruction  (NROI)  82-3,  Satellite  Debris  Mitigation-End  of  Life 
Chairman  of  the  Joint  Chiefs  of  Staff,  CJCSI  6250.01,  Satellite  Communications 

Terms 

Apogee — The  point  of  a  satellite' s  greatest  distance  from  Earth  and  minimum  velocity. 

Combatant  Command  (COCOM) — Non-transferable  command  authority  established  by  title  10, 
United  States  Code,  Section  164,  exercised  only  by  commanders  of  unified  combatant  commands. 
COCOM  is  the  authority  of  a  Combatant  Commander,  to  perform  those  functions  of  command  over 
assigned  forces  involving  organizing  and  employing  commands  and  forces,  assigning  tasks,  designating 
objectives,  and  giving  authoritative  direction  over  all  aspects  of  military  operations,  joint  training,  and 
logistics  necessary  to  accomplish  the  mission  assigned  to  the  command.  COCOM  provides  full  authority 
to  organize  and  employ  commands  and  forces  as  the  CINC  considers  necessary  to  accomplish  assigned 
missions. 

Military  Satellite  Communications  (MILSATCOM) — The  satellite  communications  resources  that  are 
owned  and  operated  by  DOD  primarily  in  the  government  frequency  bands.  USCINCSPACE  COCOM 
systems  include,  but  are  not  limited  to.  Defense  Satellite  Communications  System  (DSCS),  Fleet  Satellite 
(FLTSAT),  UHF  Follow-on  (UFO),  and  Milstar. 

Perigee — The  point  of  a  satellite’s  closest  approach  to  Earth  and  greatest  velocity. 

SATCOM  System  Expert  (SSE) — The  component  or  designated  organization  responsible  for  providing 
the  technical  planning  and  functions  in  support  of  the  operational  management  of  a  specific  SATCOM 
constellation. 

Satellite  Control  Authority  (SCA) — The  authority  to  provide  Telemetry,  Tracking  and  Commanding 
(TT&C)  of  the  satellite’s  bus  and  to  provide  control  and  management  of  the  satellite’s  payload.  Authority 
is  inherent  in  the  Operational  Control  (OPCON)  of  a  satellite  or  can  be  delegated  in  part  or  entirely  to 
another  organization. 

Secondary  Payloads — Those  payloads  placed  on  a  satellite  that  provide  a  different  capability  than  its 
primary  payload.  The  secondary  capability  of  a  satellite  system  is  taken  into  consideration  when  making 
a  decision  to  dispose  of  a  satellite, 

USSPACECOM  Components  — The  three  service  components  of  USSPACECOM  are  U.S.  Army  Space 
Command  (USARSPACECOM),  Naval  Space  Command  (NAVSPACECOM),  and  14th  Air  Force 
(SPACEAF).  14th  Air  Force,  Vandenberg  AFB,  California,  which  falls  under  the  administrative  control 
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of  Air  Force  Space  Command  (AFSPC),  Peterson  AFB,  Colorado,  is  dual  hatted  as  SPACEAF  and 
operates  as  USSPACECOM’s  functional  component. 
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Attachment  2 

METRICS  TO  MEASURE  COMPLIANCE  WITH  POLICY 

A2.1.  The  purpose  of  metrics  is  to  identify  the  established  goals  or  standards  based  on  our  customers 
needs  and  measure  how  well  we  meet  the  standards  based  on  our  customers  feedback.  The  key  metrics  in 
the  disposal  process  are  satellites  that  are  placed  within  disposal  regions  for  present  satellite  systems  and 
incorporation  of  disposal  capabilities  for  future  systems.  These  metrics  will  be  reviewed  and  revised 
every  3  years,  or  as  needed. 

A2.2.  Execute  Satellite  Placement  into  a  Disposal  Orbit.  This  metric  measures  the  number  of  satellites 
placed  into  a  proper  disposal  region  after  having  been  declared  non-mission  capable  or  subsystem  dis¬ 
posal  criteria  are  met.  The  goal  is  to  place  all  non-mission  capable  satellites  into  a  specified  disposal 
region. 

A2.2.1.  Customer:  Primary-USSPACECOM;  secondary-DOD,  National  and  commercial  users. 

A2.2.2.  Goal  and  source:  100  percent  insertion  to  proper  disposal  region,  determined  by  Cheyenne 
Mountain  Operations  Center  (CMOC). 

A2.2.3.  Data  source:  Organization  with  satellite  system  responsibility  and  CMOC. 

A2.2.4.  Calculation  method:  Percentage  determined  by  dividing  the  number  of  satellites  reaching  a 
disposal  region  by  the  total  number  of  satellites  being  disposed  of  (not  to  include  those  satellites 
undergoing  end-of-life  testing  prior  to  disposal). 

A2.3.  Satellite  Programs  Planning  for  Disposal  Capabilities.  This  metric  compares  the  number  of  new 
satellite  programs  being  developed  with  the  number  of  these  programs  that  incorporate  a  disposal  capabil¬ 
ity  into  the  system  that  meets  the  requirements  of  this  directive.  As  long  as  a  satellite  has  the  capability  to 
meet  the  disposal  criteria  upon  launch,  this  requirement  is  being  met. 

A2.3.1.  Customer:  Primary-USSPACECOM;  secondary-DOD,  National  and  commercial  users. 

A2.3.2.  Goal  and  source:  100  percent  of  future  satellite  programs  being  developed  incorporate  satel¬ 
lite  disposal  capabilities,  as  determined  by  System  Program  Offices  (SPOs). 

A2.3.3.  Data  source:  Space  and  Missile  Command  (SMC)  SPOs. 

A2.3.4.  Calculation  method:  Percentage  determined  by  dividing  the  number  of  new  programs  that 
incorporate  a  disposal  capability  by  the  number  of  new  satellite  programs. 
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